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Abstract The following chapter provides an overview of the techniques used to 
understand Self-Organised Criticality (SOC) by performing computer simulations. 
Those are of particular significance in SOC, given its very paradigm, the BTW 
(Bak-Tang-Wiesenfeld) sandpile, was introduced on the basis of a process that is 
conveniently implemented as a computer program. The chapter is divided into three 
sections: In the first section a number of key concepts are introduced, followed by 
four brief presentations of SOC models which are most commonly investigated or 
which have played an important part in the development of the field as a whole. The 
second section is concerned with the basics of scaling with particular emphasis of 
its role in numerical models of SOC, introducing a number of basic tools for data 
analysis such as binning, moment analysis and error estimation. The third section is 
devoted to numerical methods and algorithms as applied to SOC models, address- 
ing typical computational questions with the particular application of SOC in mind. 
The present chapter is rather technical, but hands-on at the same time, providing 
practical advice and even code snippets (in C) wherever possible. 



1.1 Introduction 

The concept of Self-Organised Criticality (SOC) 1 was introduced by Bak et al. 
(1987) on the basis of a computer model, the famous BTW Sandpile. The notion 
of "computer model" and "simulation" used here is subtle and can be misleading. 
Often the models are not meant to mimic a particular (natural) phenomenon, but 
are intended to capture merely what is considered to be the essential interaction ob- 
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served in a natural phenomenon. Per Bak in particular, had the tendency to name 
models according to their appearance rather than their purpose and so the "Sandpile 
Model" may not have been envisaged to display the dynamics of a sandpile. The sit- 
uation is clearer in the case of the "Forest Fire Model" (Bak et al, 1990), which was 
developed as a model of turbulence much more than as a model of fires in woods. 

In particular in the early days of SOC modelling, the models were sometimes 
referred to as "cellular automata" Olami et al. (1992); Lebowitz et al. (1990), which 
caused some consternation {e.g. Grassberger, 1994), as cellular automata normally 
have discrete states and evolve in discrete time steps according to deterministic rules 
in discrete space {i.e. a lattice). The term "coupled map lattice" (Kaneko, 1989) can 
be more appropriate for some models, such as the Olami-Feder-Christensen Model 
descusssed below (discrete space, continuous state and possibly continuous time). 

The terminology of "numerical modelling" has always been somewhat confus- 
ing. Many of the models considered in SOC do not model a natural phenomenon 
and so their numerical implementation is not a "numerical simulation" in the sense 
that they mimic the behaviour of something else. There are notable exceptions, how- 
ever, such as the Forest Fire Model (Bak et al., 1990) mentioned above and the Oslo 
ricepile model (Christensen et al., 1996). SOC models generally are not "models of 
SOC", rather they are algorithmic prescriptions or "recipes" for a (stochastic) pro- 
cess that is believed to exhibit some of the features normally observed in other SOC 
models. In that sense, the terminology of terms like "SOC models" and "simulation" 
or even "simulating an SOC model" is misleading — most of these models are not 
simplified versions or idealisations of some physical process or anything else that is 
readily identified as "SOC", but recipes to produce some of the behaviour expected 
in an SOC system. 

To this day, a large fraction of the SOC community dedicate their research to 
computer models. Initially, the motivation {e.g. Zhang, 1989; Manna, 1991) was 
to find models displaying the same universal behaviour as the BTW (Bak-Tang- 
Wiesenfeld) Sandpile. This was followed by an era of proliferation, when many 
new models, belonging to new universality classes where developed. More recently, 
in a more reductionistic spirit, new models are mostly developed to isolate the role 
of particular features and to extract and identify their effect {e.g. Tadic and Dhar, 
1997). A lot of numerical research into SOC nowadays happens "en passant", as 
SOC is identified in a model for a phenomenon that originally was not considered 
to be related to SOC {e.g. Burridge and Knopoff, 1967). 

Virtually all SOC (computer) models consist of degrees of freedom interacting 
with (nearest) neighbours located on a lattice. The degrees of freedom may be pa- 
rameterised by continuous or discrete variables, in the following denoted z n , where 
n is a position vector on the lattice. A slow, external driving mechanism (in short, 
external drive) slowly loads the system, i.e. the local variables are slowly increased, 
also referred to as "charging a site". That might happen uniformly (sometimes called 
global drive) or at individual lattice sites (sometimes called point drive). The driv- 
ing might happen at randomly chosen points or by random increments, both of 
which is in the literature referred to as random driving. The dynamics of an SOC 
model is non-linear, i.e. there is no linear equation of motion that would describe 
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their dynamics. 2 The response of the system is triggered by a local degree of free- 
dom overcoming a threshold, beyond which relaxation and thus interaction with 
other degrees of freedom and the outside world takes place. A site where that hap- 
pens is said to topple and to be active. The interaction might lead to one of the 
neighbours exceeding its threshold in turn, triggering another relaxation event. The 
totality of the relaxations constitutes an avalanche. When the avalanche has fin- 
ished, i.e. there are no active sites left, the system is in a state of quiescence. In 
SOC models, driving takes place only in the quiescent state (separation of time 
scales, below). If the external drive acts at times when an avalanche is running, it 
might lead to a continuously running avalanche (e.g. Corral and Paczuski, 1999). 

In many models the degree of freedom at every site measures a resource that 
is conserved under the dynamics. To balance the external drive, in most models 
dissipation has to take place in some form: Bulk dissipation takes place when 
the resource can get lost in the local interaction. Boundary dissipation refers to 
the situation when the resource is lost only in case a boundary site relaxes. The 
necessary flux of the resource towards the boundaries has been suggested as some 
of the key mechanisms in SOC (Paczuski and Bassler, 2000b). In some models, 
such as the Bak-Sneppen Model (Bak and Sneppen, 1993) or the Forest-Fire-Models 
(Henley, 1989; Bak et al, 1990; Drossel and Schwabl, 1992a), no (limited) resource 
can be identified and therefore the notion of dissipation and conservation is not 
meaningful. 

The question whether conservation is a necessary ingredient of SOC has driven 
the evolution of SOC models in particular during the 1990s. In fact, early theoretical 
results by Hwa and Kardar (1989a) suggested that bulk dissipation would spoil the 
SOC state. Models like the OFC Model (Olami et al, 1992, also Bak and Sneppen, 
1993; Drossel and Schwabl, 1992a) questioned that finding. Different theoretical 
views have emerged over time: Lauritsen et al.'s (1996) self-organised branching 
process (Zapperi et al, 1995) contains dissipation as a relevant parameter which 
has a limitting effect on the scaling behaviour. Juanico et al. (2007) restored the 
SOC state of the self-organised branching process by implementing a mechanism 
that compensates for the non-conservation by a "matching condition" not dissim- 
ilar from the mechanism used in the mean-field theory by Pruessner and Jensen 
(2002b). That, in turn, was labelled by Bonachela and Munoz (2009) as a form of 
tuning. More recent field-theoretic work (Pruessner, 2012b) points at conservation 
as a symmetry responsible for the cancellation of mass-generating diagrams, an ef- 
fect that may equally be achieved by other symmetries. 

The external drive, the ensuing sequence of avalanches and the evolution of the 
model from one quiescent state to the next happen on the macroscopic time scale, 
where time typically passes by one unit per avalanche. As the system size is in- 
creased, avalanches are expected to take more and more relaxations to complete. 



2 It is very instructive to ask why a non-linearity is such a crucial ingredient. Firstly, if all inter- 
actions were linear, one would expect the resulting behaviour to correspond to that of a solvable, 
"trivial" system. Secondly, linearity suggests additivity of external drive and response, so responses 
would be expected to be proportional to the drive, a rather boring behaviour, not expected to result 
in scale invariance. 
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Their duration is measured on the microscopic time scale. In the thermodynamic 
limit, i.e. at infinite system size, the infinite duration of an avalanche of the micro- 
scopic time scale and the finite driving rate on the macroscopic time scale amount 
to a complete separation of time scales. In general, the separation of time scales 
is achieved in finite systems provided that no driving takes place when any site is 
active, because the times of quiescence, measured on the microscopic time scale, 
can be thought of as arbitrarily long. As a result, the avalanching in these systems 
becomes intermittent. 

Separation of time scales is widely regarded as the crucial ingredient of SOC, 
maybe because it is conceived (and criticised as such) as a substitute of the tuning 
found in traditional critical phenomena (also Jensen, 1998). In numerical models, 
it normally enters in a rather innocent way — the system is not driven while an 
avalanche is running. This, however, requires some global supervision, a "babysit- 
ter" (Dickman et ai, 2000) or a "farmer" (Broker and Grassberger, 1999). In some 
models the separation of time scales can be implemented explicitly (Bak and Snep- 
pen, 1993) in the relaxational rule. What makes the separation of time scales very 
different from other forms of tuning is that it eliminates a dimensionful, finite scale, 
such as the frequency with which an avalanche is triggered. 3 In traditional criti- 
cal phenomena, scaling comes about due to the presence of a dimensionful, finite 
energy scale 4 , where entropic contributions to the free energy compete with those 
from the internal energy promoting order. In most SOC models, it is pretty obvi- 
ous that scaling would break down if time scales were not explicitly separated — 
avalanches start merging and eventually intermittency is no longer observed (Corral 
and Paczuski, 1999). 

SOC models are normally studied at stationarity, when all correlations origi- 
nating from the initial state (often the empty lattice) are negligible. Reaching this 
point is a process normally referred to as equilibration. The equilibration time is 
normally measured as the number of charges by the external drive required to reach 
stationarity. For some models, exact upper bounds for the equilibration time are 
known (Dhar et ai, 1995; Corral, 2004a; Dhar, 2004, e.g.). In deterministic models, 
a clear distinction exists between transient and recurrent states, where the former 
can appear at most once, and the latter with a finite frequency provided the number 
of states overall is finite. In fact, this frequency is the same for all recurrent states, 
depending on the driving, which can be at one site only or randomly and indepen- 
dently throughout. A detailed proof of such properties can be cumbersome (Dhar, 
1999a,b). 

The statistics of the avalanches, their size as well as their extent in space and 
in time, is collected and analysed. SOC is usually said to be found in these mod- 
els when the statistics displays a scaling symmetry, governed by only one upper 
cutoff which diverges with the system size. In principle, a Gaussian possesses this 



3 In the field theory of SOC, the cancellation of diagrams occurs precisely when stationarity is 
imposed for the density of particles resting (and their correlations) in the limit (O — » 0, i.e. in the 
long time limit. 

4 For example kgT c in the Ising Model (Stanley, 1971). 
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scaling symmetry, 5 but not a single important SOC model has a Gaussian event size 
distribution. On the contrary, the avalanche statistics of all models discussed be- 
low deviates dramatically from a Gaussian, thus suggesting that avalanches are not 
the result of essentially independent patches of avalanching sites creating a bigger 
overall avalanche. Rather, sites are strongly interacting, thereby creating the over- 
all event. The purpose of numerical simulations is to characterise and quantify this 
interaction and its effect, as well as extracting universal quantities, which can be 
compared with those found in other systems. 



1.1.1 Observables 

As for the methods of analysis, they have matured considerably over the past 
decades. The initial hunt for 1 // noise in temporal signals has given way to the study 
of event size distributions. As a matter of numerical convenience, these distributions 
are often characterised using moments, some of which are known exactly. Since 
the beginning of computational physics, moments and cumulants have been the 
commonly used method of choice to characterise critical phenomena (Binder and 
Heermann, 1997). It is probably owed to the time of the late 1980's that memory- 
intensive observables such as entire distributions became computationally afford- 
able and subsequently the centre of attention in SOC. 

To this day, the analysis of moments in SOC is still often regarded as an unfor- 
tunate necessity to characterise distributions, which are difficult to describe quan- 
titatively. Apart from the historic explanation alluded to above, there is another, 
physical reason for that, the avalanche size exponent T. In traditional critical phe- 
nomena, the corresponding exponent of the order parameter distribution is fixed at 
unity in the presence of the Rushbrooke and the Josephson scaling law (Christensen 
et at, 2008). The deviation of t from unity, which implies that the expected event 
size does not scale like the characteristic event size, is another distinctive feature 
of SOC. To some extent, the exponent T can be extracted from the avalanche size 
distribution (almost) by inspection. In a moment analysis, on the other hand, it is 
somewhat "hidden" in the details. 

The most important observables usually extracted from an SOC model are thus 
the scaling exponents, such as x, D (avalanche dimension), a (avalanche duration 
exponent) and z (dynamical exponent) discussed below. Here, the two exponents 
D and z are generally regarded as more universal than T and a, as the former is of- 
ten "enslaved" by an exact scaling law related to the average avalanche size, and the 
latter by a similar scaling law based on the "narrow joint distribution assumption", 
discussed in Sec. 1.2. Generally, all observables that are universal or suspected to 
be are of interest. This includes the scaling function (Sec. 1 .2) which is most easily 
characterised by moment ratios, corresponding to universal amplitude ratios, tradi- 

5 The basic example &(s) = s~ 1 @(s/s c ) with Sf(jc) = 2*exp(— x 2 )/^Jk is normalised and has 
avalanche size exponent z = 1, as defined in Eq. (1.3). 'Without the pre-factor x in Sf(jc) the graph 
looks surprisingly similar to a PDF as typically found in SOC models. 
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tionally studied in equilibrium critical phenomena (Privman et ai, 1991; Salas and 
Sokal, 2000). 



1.1.2 Models 

There is wide consensus on a number of general features of SOC models which seem 
to play a role in determining the universality class each belongs to. The very first 
SOC model, the BTW model, was essentially deterministic, i.e. there was no ran- 
domness in the bulk relaxation. A given configuration plus the site being charged 
next determines the resulting configuration uniquely. Even in these models, how- 
ever, there can be a degree of stochasticity, namely when the site to be charged by 
the external drive is chosen at random. Finally, even when this is not the case, i.e. 
external drive and internal relaxation are deterministic, initial conditions are often 
chosen at random and averaged over. 

Deterministic SOC models have the great appeal that they are "autonomous" (in 
a non-technical sense) or "self-sufficient" in that they do not require an additional 
source of (uncorrelated) noise. It is difficult to justify the existence of an external 
source which produces white, Gaussian noise, as that noise correlator, (j](t)r](t')') = 
2r 2 8(t — t'), itself displays a form of scaling (rj(at)rj(at')) = a~ l (~n{t)ri(t')). 
The presence of an external (scaling) noise source seems to demote an SOC model 
to a conversion mechanism of scale invariance, which becomes most apparent when 
the respective model is cast in the language of stochastic equations of motion, i.e. 
Langevin equations. 

Famous examples of deterministic SOC models, which do not require an external 
noise source for the relaxation process, are the BTW model with deterministic drive 
(Bak et ai, 1987, but Creutz, 2004), the OFC model (Olami et ai, 1992) and, closely 
related, the train model (de Sousa Vieira, 1992). Of these only the latter has been 
studied extensively in the absence of all stochasticity. 

Most SOC models, however, have a strong stochastic component, i.e. there is 
some randomness in the relaxation mechanism that gives rise to avalanches. In fact, 
models with some form of built-in randomness seem to give cleaner scaling be- 
haviour, suggesting that deterministic models get "stuck" on some trajectory on 
phase space, where some conservation law prevents them from exploring the rest 
of phase space (Bagnoli et ai, 2003; Casartelli et ai, 2006). Notably, randomis- 
ing the BTW model seems to push it into the Manna universality class (Karmakar 
et ah, 2005). The latter model is probably the simplest SOC model displaying the 
most robust and universal scaling behaviour (Huynh et ai, 201 1). Due to the noise, 
trajectories of particles deposited by the external drive are those of random walkers. 

The second dividing line distinguishes Abelian and non-Abelian models. The 
term was coined by Dhar (1990) introducing, strictly speaking, the Abelian Sandpile 
Model, by re-expressing the original BTW Model (Bak et ai, 1987) in terms of units 
of slope rather than local particle numbers. This convenient choice of driving and 
boundary conditions renders the model unphysical as entire rows of particles are 
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added and removed at once. At the same time, however, the model's final state after 
two consecutive charges at two different sites becomes independent from the order 
in which the charges and the subsequent relaxations are carried out. Practically all 
analytical insight into the BTW model is based on Dhar's (1990) Abelian version. 
Because it is easier to implement, it has also favoured in numerical simulations. 

The term "Abelian" seems to suggest the existence of a (commutative) group, i.e. 
a set of operators closed under consecutive application, associative and containing 
inverse and an identity. For most SOC models referred to as Abelian, no such group 
is known, for example because operators do not exist explicitly, or the associative 
property makes little sense, similarly for the identity. Crucially, inverse operators 
rarely exist. To label a model Abelian therefore normally means that the final state 
does not depend on the order in which external charges are applied, i.e. the model 
updating operators (whether or not they exist), which drive it at various locations, 
commute. Because the final state is unique only in the case of deterministic mod- 
els, stochastic models are Abelian provided that the statistics of the final state does 
not depend on the order in which external charges are applied (Dhar, 1999b). The 
operators, which generally depend on the site the driving is applied to, of deter- 
ministic models apply to a model's state and take it from one quiescent state to the 
next. The operators in a stochastic model act on the distribution of states, i.e. they 
are the Markov operators. A deterministic model can be cast in the same language, 
however, the Markov operators then correspond to simple permutation matrices. 

While Abelianness originally refers to the evolution of a model on the macro- 
scopic time scale, it is generally used to characterise its behaviour on the micro- 
scopic timescale, i.e. the step-by-step, toppling-to-toppling update. It is therefore 
usually concluded that the properties of avalanches and their statistics is indepen- 
dent from the order of toppling of multiple active sites. 

Strictly, however, the Abelian symmetry does not apply to the microscopic time 
scale, at least for two reasons. Firstly, the Abelian operators apply, a priori, only to 
the avalanche-to-avalanche evolution, i.e. the macroscopic time scale. What is more, 
they apply to the final state and its statistics, but not necessarily to the observables. 
Applying charges at two different sites of an Abelian SOC model, starting from 
the same configuration, results in the same final state (or its statistics) regardless 
of the order in which the charges were applied, but not necessarily in the same 
pair of avalanche sizes produced. On the basis of the proof of Abelianness, at least 
in deterministic models, this limitation is alleviated by the insight that the sum of 
the avalanche sizes is invariant under a change of the order in which the model is 
charged. 

As for the second reason, many models come with a detailed prescription of the 
microscopic updating procedure and therefore the microscopic time scale. Strictly, 
the invariance under a change of order of updates on the microscopic time scale 
thus applies to different models. The situation corresponds to equating different 
dynamics in the Ising model: For some observables, Glauber dynamics is different 
from Heat Bath dynamics, yet both certainly produce the same critical behaviour. In 
fact, choosing different dynamics (and thereby possibly introducing new conserved 
symmetries) can lead to different dynamical critical behaviour. 
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Revisting the proof of Abelianness, however, generally reveals that the caveats 
above are overcautious. The very proof of Abelianness on the macroscopic time 
scale uses and develops a notion of Abelianness on the microscopic time scale. 
This connection can be made more formally, once it has been established that any 
configuration, quiescent or not, can be expressed by applying a suitable number of 
external charges on each site of an empty lattice. 

Abelianness generally plays a major role in the analytical treatment of SOC mod- 
els, because it allows significant algebraic simplifications, not least when the dynam- 
ics of a model is written in terms of Markov matrices. It applies, generally, equally 
to recurrent and transient states, where no inverse exists. It remains highly desirable 
to demonstrate Abelianness on the basis of the algebra, once that is established as a 
suitable representation of a model's dynamics. 

In the following section a few paradigmatic models of SOC are introduced: The 
BTW Model, the Manna Model, the OFC Model and the Forest Fire Model. 



1.1.2.1 The BTW Model 

The BTW Model was introduced together with the very concept of SOC (Bak et ah, 
1987), initially to explain the "ubiquity" of l/f noise. Of course, since then, SOC 
has been studied very much in its own right. Like virtually all SOC models, the 
BTW Model consists of a set of rules that prescribe how a local degree of freedom 
Zi on a (/-dimensional lattice with sites ; is to be updated. There are two different 
stages, namely the relaxation and the driving, the latter considered to be slow com- 
pared to the relaxation, i.e. the relaxation generally is instantaneous and never occurs 
simultaneously with the driving (separation of time scales). In the Abelian version 
of the BTW Model (Dhar, 1990), the driving consists of adding a single slope unit 
(Kadanoff et ah, 1989) to a site, that is normally picked uniformly and at random. 
The lattice is often initialised with zi = for all i. 

If the driving leads to any of the Zi exceeding the the critical slope z c (also referred 
to as the critical height or threshold, depending on the view) at a site / a toppling 
occurs whereby Zi is reduced by the coordination number q of the site and Zj of 
every nearest neighbour j increases by one (sometimes referred to as charging). In 
principle both q and z c can vary from site to site and such generalisations are trivial 
to implement. It is common to choose z c = q — 1 . 

The rules of the BTW Model can be summarised as follows: 

Initialisation: All sites i are empty, zi = 0. 

Driving: One unit is added at a randomly chosen (or sometimes fixed) site i, i.e. 

Zi—Zi + l. 

Toppling: A site with z, > z c = q — 1 (called active) distributes one unit to the q 
nearest neighbouring sites j, so that zi —* Zi~q and zj — ► Zj+ 1. 

Dissipation: Units are lost at boundaries, where toppling site ; loses q units, zi — » 
Zi — q, yet less than q nearest neighbours exist, which receive a unit. 

Time progression: Time progresses by one unit per parallel update, when all ac- 
tive sites are updated at once. 
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A toppling can trigger an avalanche, as charged neighbours might exceed the 
threshold in turn, possibly by more than one unit. Strictly, the BTW Model is up- 
dated in parallel, all sites topple at once whose local degree of freedom exceeds 
the threshold at the beginning of a time step. Microscopic time then advances by 
one unit. This way, zt might increase far beyond z° before toppling itself. As long 
as zi > z c for any site i, the sites in the model carry on toppling. The totality of 
the toppling events is an avalanche. In the Abelian BTW model as refined by Dhar 
(1990), the final state of the model does not depend on the order in which external 
charges are applied. In the process of the proof of this property, it turns out that 
the order of processing any charges during the course of an avalanche neither af- 
fects the final state nor the size of the avalanche triggered. Using a parallel updating 
scheme or not therefore does not change the avalanche sizes recorded. As the order 
of updates defines the microscopic time scale, a change in the updating procedure, 
however, affects all observables dependent on that time, such as avalanche duration 
or correlations on the fast time scale. 

To keep the prescription above consistent with the notion of boundary sites, 
where toppling particles are to be lost to the outside, boundary sites have to be 
thought of as having the same number of nearest neighbours as any other, equiva- 
lent site in the bulk, except that some of their neighbours are not capable of toppling 
themselves. For numerical purposes it is often advisable to embed a lattice in some 
"padding" (a neighbourhood's "halo", see Sec. 1.3.2.2, p. 39), i.e. sites that cannot 
topple but are otherwise identical to all other sites. 

The sum of the slope units residing on a given site i and those residing on its 
nearest neighbours remains unchanged by the toppling of site i, i.e. the bulk dynam- 
ics in the BTW are conservative. Dissipation occurs exclusively at the boundary and 
every slope unit added to the system in the bulk must be transported to the boundary 
in order to leave the system. 

The original version of the BTW model is defined in terms of local heights, so 
that the height differences give rise to the slope zi, which has to reach q in order to 
trigger an a toppling. While this is a perfectly isomorphic view of the BTW, driving 
it in terms of height units has a number of unwanted implications. In particular, 
it loses its Abelianness. For that reason, the original version of the BTW is rarely 
studied numerically nowadays. 

The BTW Model is deterministic apart from the driving, which can be made 
deterministic as well, simply by fixing the site that receives the external charge that 
triggers the next avalanche. Even when slope units do not move independently at 
toppling, a randomly chosen slope unit being transported through a BTW system 
describes the trajectory of a random walker trajectories (Dhar, 1990), essentially 
because every possible path is being realised (just not independently, but all with the 
correct weight). As a result, the average avalanche size (s) can be calculated exactly; 
The number of moves a slope unit makes on average from the time of being added 
by the external drive to the time it leaves the system through an open boundary is 
equal to the expected number of charges it causes. The expected number of charges 
(caused by the movement of all slope units taking part in an avalanche) per slope 
unit added is thus exactly equal to the expected number of moves a slope unit makes 
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until it leaves the system, i.e. its escape time. If the avalanche size is measured by 
the number of topplings, which is more common, the expected number of moves 
has to be divided by the number of moves per toppling, q in the present case. Higher 
moments of the avalanche size, or, say, the avalanche size conditional to non-zero 
size (i.e. at least one site toppling in every avalanche), cannot be determined using 
the random walker approach, as they are crucially dependent on the interaction of 
toppling sites. 

Due to the random walker property of the slope units added, the scaling of the 
average avalanche size thus merely depends on the particularities of the driving. If 
the driving is random and uniform, then (s) ccL 2 for any of-dimensional hypercubic 
lattice and (Ruelle and Sen, 1992) 

(s)=^(L+l)(L + 2) (1.1) 

in one dimension with two open boundaries, where the avalanche size is the number 
of topplings per particle added. However, the dynamics of the BTW Model in one 
dimension is trivial, so that the model is usually studied only in d = 2 and beyond. 

Because (or despite of) its deterministic nature, a large number of analytical re- 
sults are known, in one dimension (Ruelle and Sen, 1992) but more importantly in 
two dimensions (Majumdar and Dhar, 1992), not least on the basis of (logarithmic) 
conformal field theory (e.g. Majumdar and Dhar, 1992; Ivashkevich, 1994; Mahieu 
and Ruelle, 2001; Ruelle, 2002; Jeng, 2005). Unfortunately, to this day, the scaling 
of the avalanche size distribution in dimensions d > 2 remains somewhat unclear. 
Numerically, results are inconclusive, as different authors quote widely varying re- 
sults for d = 2 (Vespignani and Zapperi, 1995; Chessa et al., 1999a; Lin and Hu, 
2002; Bonachela, 2008, e.g.), possibly due to logarithmic corrections (Manna, 1990; 
Liibeck and Usadel, 1997; Liibeck, 2000) 

A major insight into the collective dynamics of toppling sites was the decompo- 
sition of avalanches into waves (Ivashkevich et al, 1994), which was later used by 
Priezzhev et al. (1996) to conjecture t = 6/5 for the avalanche size exponent in two 
dimensions. No site in an avalanche can topple more often than the site at which 
the avalanche was triggered. Not allowing that first site to topple therefore stops 
the avalanche from progressing any further and each toppling of the first site thus 
defines a wave of toppling. 

While the BTW Model has been crucial for the formation of the field of SOC as 
a whole, its poor convergence beyond one dimension has made it fall in popularity. 
One may argue that the determinism of the dynamics is to blame, as found in other 
models (Middleton and Tang, 1995). Indeed, adding some stochasticity makes the 
BTW Model display the universal behaviour of the Manna Model discussed in the 
next section (Cernak, 2002; Cernak, 2006). 

The exponents reported for the BTW Model vary greatly. In two dimensions, the 
value of T found in various studies ranges from 1 (Bak et al., 1987) to 1.367 (Lin 
and Hu, 2002) and that for D from 2.50(5) (De Menech et al., 1998) to 2.73(2) 
(Chessa et al., 1999a). Similarly a is reported from 1.16(3) (Bonachela, 2008) to 
1.480(11) (Liibeck and Usadel, 1997) and z from 1.02(5) (De Menech and Stella, 
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2000) to 1.52(2) (Chessa et ai, 1999a). Using comparatively large system sizes, 
Dorn et al. (2001) found exponents that seem to vary systematically with the system 
size with little or no chance to identify an asymptotic value. 

The first exactly solved SOC model was the Dhar-Ramaswamy Model (Dhar and 
Ramaswamy, 1989) which is the directed variant of the BTW Model. The direct- 
edness means that during an individual avalanche, sites are never re-visted, which 
effectively suppresses spatial correlations. Random drive of the model results in a 
product state, where sites taking part in an avalanche form a "compact" patch {i.e. 
they have no holes), which is delimited by boundaries describing a random walk. 
The exponents in d = d± + 1 dimensions are given analytically by D = 1 + d±/2, 
D(2-t) = l,z = 1 and£>(T- 1) = z(a- 1), which implies a = D and x = 2- 1/D 
(Dhar and Ramaswamy, 1989; Christensen, 1992; Christensen and Olami, 1993; 
Tadic and Dhar, 1997; Kloster et al., 2001). For example, in d = 1 + 1 dimensions 
(directed square lattice), exponents are D = 3/2, x = 4/3, z = 1 and a = 3/2. Mean- 
field exponents apply at d = 2 + 1 and above. 

1.1.2.2 The Manna Model 

The Manna (1991) Model was originally intended as a simplified version of the 
BTW Model but has since then acquired the status of the paradigmatic representative 
of the largest (and maybe the only) universality class in SOC, generally refered to as 
the Manna, Oslo (Christensen et ai, 1996) or C-DP (conserved directed percolation, 
Rossi et ai, 2000) universality class. 

The Manna Model displays robust, clean critical behaviour in any dimension 
d 5= 1, characterised by non-trivial exponents below d = 4 (Liibeck and Heger, 
2003b). Originally, it is defined as follows: The external drive adds particles at ran- 
dom chosen sites i, i.e. the local degree of freedom increases by one, z, — ► zi+ 1. 
If a site exceeds the threshold of z c = 1 it topples, so that all its particles are redis- 
tributed to the nearest neighbours, which are chosen independently at random. After 
the toppling of site i, the local degree of freedom is therefore set to z, = 0, while 
the total increase of the zj at the nearest neighbours j of i maintains conservation. 
Again, as in the BTW model, non-conservation at boundary sites can be thought of 
as been implemented by sites that never topple themselves. 

Charging neighbours might push their local degree of freedom beyond the thresh- 
old and they might therefore topple in turn. When a site topples, all particles present 
there at the time of toppling are transferred to its neighbour (maybe to a single one) 
and it is therefore crucial to maintain the order of (parallel) updates. The model is 
thus non-Abelian. In fact, the notion of Abelianness was initially restricted to de- 
terministic models (Milshtein et ai, 1998). However, Dhar (1999a) introduced a 
version of the Manna Model which is Abelian in the sense that the statistics of the 
final state remains unchanged if two consecutive external charges (by the driving) 
are carried out in reverse order. In that version of the Manna Model, a toppling site 
redistributes only 2 of its particles, i.e. the number of particles redistributed at a top- 
pling does not depend on z, itself. The difference between the BTW Model and the 
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Manna Model lies thus merely in the fact that only two particles are re-distributed 
when a site topples in the Manna Model (irrespective of the coordination number of 
the site) and that the receiving sites are are picked at random. 
In summary, the rules of the Abelian Manna Model are: 

Initialisation: All sites i are empty, zi = 0. 

Driving: One unit is added at a randomly chosen (or sometimes fixed) site i, i.e. 
Zi -^Zi + 1. 

Toppling: A site with z, > z c = 1 (called active) distributes one unit to 2 randomly 
and independently chosen nearest neighbouring sites j, so that z, — » z, ■ — 2 and 

Zj—'Zj+l. 

Dissipation: Units are lost at boundaries, where the randomly chosen nearest 

neighbour might be outside the system. 
Time progression: Originally, time progresses by one unit per parallel update, 

when all active sites are updated at once. 

That the scaling in one dimension is not as clean as in higher dimension may be 
caused by logarithmic corrections (Dickman and Campelo, 2003). Nevertheless, it 
has been possible to extract consistent estimates for exponents in dimensions d = 1 
to d = 5 (Liibeck and Heger, 2003b; Huynh et al, 2011; Huynh and Pruessner, 
2012). Because some of its exponents are so similar to that of the directed perco- 
lation universality class (Janssen, 1981; Grassberger, 1982; Hinrichsen, 2000) there 
remains some doubt whether the Manna Model really represents a universality class 
in its own right (Munoz et al, 1999; Dickman et al, 2002). The problem is more 
pressing in the fixed energy version (Dickman et al, 1998; Vespignani et al, 1998) 
of the Manna Model (Basu et al, 2012), where dissipation at boundaries is switched 
off by closing them periodically, thereby studying the model at a fixed amount of 
particles. The term "fixed energy sandpile" was coined to stress the conserved nature 
of the relavent degree of freedom (which may be called "energy") and to suggest a 
similar distinction as in the change of ensemble from canonical to microcanonical. 
Bonachela and Munoz (2007) suggested to study the model with different boundary 
conditions which have an impact on the Manna Model that is distinctly different 
from that on models in the directed percolation universality class. 

Because of its fixed energy version, the Manna Model is frequently studied for 
its links to absorbing state (AS) phase transitions (Dickman et al, 1998; Vespignani 
et al, 1998; Hinrichsen, 2000; Henkel et al, 2008). In fact, it has been suggested 
that SOC is due to the self-organisation to the critical point of such an AS phase 
transition (Tang and Bak, 1988; Dickman et al, 1998; Vespignani et al, 1998), 
whereby strong activity leads to a reduction of particles by dissipation, making the 
system in-active, while quiescence leads to activity due to the external drive. One 
may argue that such a linear mechanism cannot produce the desired universal critical 
behaviour without finely tuning the relevant parameters (Pruessner and Peters, 2006, 
2008; Alava et al, 2008). 

A number of theoretical results are available for the Manna Model (Vespignani 
et al, 1998, 2000; Rossi et al, 2000; van Wijland, 2002; Ramasco et al, 2004), yet 
an e-expansion (Le Doussal et al, 2002) for the Manna universality class is avail- 
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able only via the mapping (Paczuski and Boettcher, 1996; Pruessner, 2003) of the 
Oslo Model (Christensen et ai, 1996), which is the same universality class (Nakan- 
ishi and Sneppen, 1997) as the Manna Model, to the quenched Edwards- Wilkinson 
equation (Bruinsma and Aeppli, 1984; Koplik and Levine, 1985; Nattermann et ah, 
1992; Leschhorn et ai, 1997). Quenched noise and disorder are, however, notori- 
ously difficult to handle analytically. It is thus highly desirable to develop a better 
theoretical understanding of the Manna Model in its own right, including its mech- 
anism of self-organisation, and to derive an e-expansion for its exponents. 

Although the Manna Model is more frequently studied in one dimension, for 
comparison with the BTW Model above, the exponents listed in the following were 
determined numerically in two dimensions for the Abelian and the non-Abelian 
(original) variant of the Manna Model. For x they range from 1 .25(2) (Biham et ai, 
2001)to 1.28(2) (Manna, 1991; Liibeck and Heger, 2003a), for D from 2.54 (Ben- 
Hur and Biham, 1996) to 2.764(10) (Liibeck, 2000), for a from 1.47(10) (Manna, 
1991) to 1.50(3) (Chessa et ai, 1999b; Liibeck and Heger, 2003a) and for z from 
1.49 (Ben-Hur and Biham, 1996) to 1.57(4) (Alava and Munoz, 2002; Dickman 
et ai, 2002), generally much more consistent than in the BTW Model. 

As in the BTW Model, various directed variants of the Manna Model which 
are exactly solvable for similar reasons as in the deterministic case have been ex- 
tensively studied (Pastor-Satorras and Vespignani, 2000b, a; Hughes and Paczuski, 
2002; Pan et ai, 2005; Jo and Ha, 2008). They have been characterised in detail by 
Paczuski and Bassler (2000b) and related to the deterministic directed models by 
Bunzarova (2010). Exponents generally follow D = 3/2 + d±/4, which can be in- 
terpreted as the diffusive exploration of a random environment. Again, correlations 
are suppressed as sites are never re-visited in the same avalanche. As in the deter- 
ministic case, z = 1 and D(2 — t) = 1 and D(x — 1) = z(a — 1) result in D = a. In 
d = 1 + 1 exponents are x = 10/7, D = 7/4, a = 7/4 and z = 1. 

1.1.2.3 The Forest Fire Model 

The Forest Fire Model has an interesting, slightly convoluted history. Two distinct 
versions exist, which share the crucial feature that the bulk dynamics is not conser- 
vative. In the original version introduced by Bak etal. (1990) sites ;, most frequently 
organised in a (two-dimensional) square lattice with periodic boundary conditions, 
can be in one of three states Oj e {T,F,A}, corresponding to occupation by a Tree, 
by Fire or by Ash. As time t advances in discrete steps, the state changes cyclically 
under certain conditions: A Tree turns into Fire at time f + 1 if a nearest neighbour- 
ing site was on Fire at time t . In turn, a Fire at time t becomes Ash in time t + 1, and 
a site covered in Ash at time t might become occupied by a Tree at time t + 1 due to 
a repeated Bernoulli trial with (small) probability p. Starting from a lattice covered 
in trees, a single site is set on fire and the system evolves under the rules described. 
The key observable is the number of sites on fire as a function of time. 

Initialisation: All (many) sites i contain a tree (otherwise ash), c,- = T, and (at 
least) one site is on fire, a, = F . 
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Fig. 1.1 Realisation of the original Forest Fire Model by Bak et al. (1990). Ash is marked by a 
white site, Trees are black and Fires grey. 



Driving: With (small) probability p, a site i containing ash at the beginning of 
time step t contains a tree, a,- = A —* T at time t + 1. 

Toppling: A site i that contains a tree at begining of time step t and has at least one 
nearest neighbour on fire, turns into fire as well, a, = T — » F. Simultaneously, a 
site on fire at t turns into ash, a, ■ = F — > A. 

Dissipation: trees grow slowly in Bernoulli trials and are removed in the "top- 
pling". Their number is not conserved under any of the updating. 

Time progression: Time progresses by one unit per parallel update. 

The original Forest Fire Model (FFM) just described possesses an absorbing 
state from which it cannot recover within the rules given. If the fire stops spreading 
because the last site on fire is surrounded by ash, the only transition that can and 
will take place eventually occupies every site by a tree. Bak et al. (1990) originally 
suggested that occasional re-lightning might be necessary — in fact, if p is large 
enough, on sufficiently large lattices, there will always be tree to burn available. 
This, however, points to a fundamental shortcoming, as quantified by Grassberger 
and Kantz (1991), namely that the lengthscale of the relevant features of the FFM are 
determined by p. Typically, at small p, some large spiral(s) of fire keeps sweeping 
across the lattice. If p is chosen too small, the spatial extent of the spiral becomes too 
large compared to the size of the lattice and the fire eventually goes out. However, 
if a control parameter determines the characteristic length scale of the phenomenon, 
it cannot be bona fide SOC (e.g. Bonachela and Munoz, 2009). Figure 1 . 1 shows an 
example of the structures, most noticeable the fire fronts, developing. 

The name "Forest Fire Model" should be taken as a witty aide-memoire. Bak 
et al. (1990) designed the model to understand scale free dissipation with uniform 
driving as observed in turbulent flow. The model should therefore be considered 
much more as a model of turbulence that happened to look like fires spreading in 
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a forest. In the present model, perpetual fires spread across trees as they re-grow, 
which is a rather unrealistic picture; most fires in real forests are shaped by fire 
brigades, geographical and geological features and other environmental character- 
istics, as well as policies. Nevertheless, the original FFM as well as the version by 
Drossel and Schwabl (1992a), attracted significant attention as an actual model of 
forest fires, as well as other natural and sociological phenomena (Turcotte, 1999). 

There are two distinguishing features that set the FFM apart from many other 
SOC models. Firstly, the separation of time scales is incomplete, because driving 
the system by supplying new trees is a process running in parallel to the burning 
as fire spreads. Although the time scale of tree growth, parameterised by p, can in 
principle be made arbitrarily slow, the fire has to be constantly fed by new trees and 
cannot be allowed to go out, because there is no explicit re-lighting. In other words, 
the tree growth rates that still sustain fire are bounded from below. As a result, there 
are no distinct avalanches, as found in the BTW and the Manna Models. 

More importantly, however, the FFM is different from other models because it 
is non-conservative at a fundamental level. No quantity is being transported to the 
boundaries and the local degree of freedom changes without any conservation. 6 At 
the time of the introduction of the FFM, it challenged Hwa and Kardar's (1989a) 
suggested mechanism of SOC that relied on a conservation law to explain the ab- 
sence of a field-theoretic mass in the propagator. 

Other dissipative models, like the SOC version of the "Game of Life" (Bak et ah, 
1989a), the OFC model discussed in the next section (Olami et al., 1992) and the 
Bak-Sneppen Model (Bak and Sneppen, 1993) chipped away from the conservation 
argument put forward by Hwa and Kardar (1989a, 1992), Grinstein et al. (1990) and 
Socolar et al. (1993). The latter seem to have been caught by surprise by the advent 
of a variant of the FFM by Drossel and Schwabl (1992a) discussed in the following. 

The Drossel-Schwabl Forest Fire Model (DS-FFM), as it is now normally re- 
ferred to, was originally introduced by Henley (1989). It changes the original Forest 
Fire Model in two very important points: Firstly, the separation of time scales be- 
tween burning and growing is completed, so that patches of (nearest neighbouring) 
trees are burned down instantly compared to all other processes. Because fires there- 
fore burn down completely before anything else can happen, fires are set, secondly, 
explicitly by random, independent uniform lightning. The key-observables of the 
DS-FFM are the geometrical features of the clusters burned down, such as the num- 
ber of occupied sites (the mass) and the radius of gyration. 

While trees grow with rate p on every empty site {i.e. one containing ash), light- 
ning strikes with much lower rate / on every site. If it contains a tree, the fire eradi- 
cates the entire cluster of trees connected to it by nearest neighbour interactions. In 
summary: 



6 It is difficult to make the statement about non-conservation strict. After all, the state of each 
site is meant to change and allowing for that, it is always possible to trace the appearance and the 
disappearance of something back to some influxes and outfluxes. Here is an attempt in the present 
case: While the increase in the number of trees can be thought of as being due to a corresponding 
influx, they can disappear with an enormous rate by spreading fire without explicit outflux on that 
timescale. 
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Fig. 1.2 Realisation of the Drossel-Schwabl Forest Fire Model (Drossel and Schwabl, 1992a). 
Ash is marked by a white site, Trees are black. 



Initialisation: All sites i contain ash, c, = A. 

Driving: With (small) probability p, a site i containing ash at the beginning of 
time step t contains a tree, c,- = A — » T at time t+1. 

Toppling: With probability / « p, a site containing a tree at the beginning of 
time step t and the entire cluster of trees connected to it by nearest neighbour 
interactions is changed to ash, <7,- = T — > A. 

Dissipation: trees grow slowly in Bernoulli trials and are removed in the "top- 
pling". They are not conserved in any of the updates. 

Time progression: Time progresses by one unit per parallel update, toppling is 
instantaneous relative to growing trees. 

As a result entire patches of forest disappear at a time, which are re-forested 
with the same Poissonian density p. This process results in a patchy structure with 
individual islands having roughly homogeneous tree-density, Figure 1.2. 

In a change of perspective, the processes parameterised by p and / are tree 
growth attempts and lightning attempts which fail if the site is already occupied 
by a tree or does not contain one, respectively. The original definition by Drossel 
and Schwabl (1992a) still used discrete time, so that both p and / were probabil- 
ities, rather than Poissonian rates, which can be recovered by rescaling p and / 
simultaneously. However, it is common (e.g. Clar et al., 1996) to rescale time so 
that p = 1 (enforced growth on randomly picked empty sites) and to attempt p/f 
times to grow a tree before attempting to set one alight. In order to see scale-free 
cluster size distributions, a second separation of timescales is needed, whereby the 
ratio p/f diverges. 

Many of the properties of the DS-FFM are percolation-like. If it were not for 
the correlations in the tree-density, which develop because of "synchronous, patchy 
re-forestation", i.e. if the tree-density was homogeneous, then the DS-FFM would 
be a form of percolation. In particular, the cluster size distribution (of the patches 
removed and the totality of all patches present) was given by that of (well-known) 
static percolation. 
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The DS-FFM does not suffer from the same short-coming as the original FFM 
of having a well-understood typical (spiral) structure, whose size is determined by 
the single control parameter p, yet it still has one control parameter which needs to 
be finely tuned in accordance with the system size. This parameter is p/f — if it 
is too large, then the lattice will be densely filled with trees before lightning strikes 
and removes almost all of them, leaving behind essentially a clean sheet with a few 
remaining (small) islands of trees. If p/f is too small, then no dense forest ever 
comes into existence and the cluster size distribution has a cutoff not determined by 
the system size, but by that parameter. 

In extensive numerical studies (Grassberger, 2002; Pruessner and Jensen, 2002a, 
2004), the system sizes were chosen big enough for each p/f that finite size ef- 
fects were not visible, i.e. for each p/f convergence of the cluster size distribution 
&(s;L) in the system size L was achieved. However, these studies revealed that the 
DS-FFM does not display simple scaling in s c = s c (p/f), Eq. (1.3) (Sec. 1.2.1). 
While £P(s) /s~ T converges in the thermodynamic limit (as it should, trivially) for 
any T, there is no choice of T so that the remaining functional profile depends only on 
the ratio s/s c (p/f). Instead, £?(s) /s~ T depends explicitly on both s and s c (p/ f), or, 
for that matter, p/f. The only feature that may display some convergence (Pruessner 
and Jensen, 2002a) is the bump in the probability density function (PDF) towards 
large s. For some choice of T, there is a small region, say [s c (p/f)/2,s c (p/f)], 
where 3^{s) /s~ x traces out a very similar graph, as if the lower cutoff sq itself was 
a divergent multiple of the upper cutoff. 7 

One may hope that finite size scaling can be recovered, taking the limit of large 
p/f and considering £P(s) /s~ T as a function of L. However, it is clear that the PDF 
trivialises in this limit, 



as the lattice is completely covered in trees before they all get completely removed 
in a singly lightning. 

Interestingly, the lack of scaling in finite s c (p/f) is not visible in the scaling 
of the moments <(*") because they are sensitive to large event sizes (at any fixed 
n > X — 1), rather than the smaller ones around the lower cutoff, whose divergence 
violates simple scaling. 

As in the BTW Model, exponents reported for the DS-FFM (if they are reported 
at all) display a fairly wide spread. In two dimensions, they are t from 1 (Drossel 
and Schwabl, 1992a) to 1.48 (Patzlaff and Trimper, 1994) and D from 1 (Drossel 
and Schwabl, 1992a) to 1.17(2) (Henley, 1993; Honeckerand Peschel, 1997). 



7 If s c (p/f) marks roughly the maximum of the bump, the PDF drops off beyond it so quickly, that 
next to nothing is known of 8?(s) beyond s c . In principle, however, if there is approximate coinci- 
dence on [s c (p/f)/2,s c (p/f)], there should also be approximate coincidence on [s c (p/f)/2,oo). 




(1.2) 
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1.1.2.4 The OFC Model 

To this day, the Olami-Feder-Christensen Model (OFC Model Olami et al, 1992) 
is one of the most popular and spectacular models of SOC. It is a simplified ver- 
sion of the Burridge-Knopoff Model (Burridge and Knopoff, 1967) of earthquakes, 
it has a tunable degree of non-conservation (including a conservative limit) with a 
clear physical meaning, it has been extensively analysed, both in time and space, 
for the effect of different boundary conditions (Middleton and Tang, 1995), and its 
one-dimensional variant (de Sousa Vieira, 1992) has been linked to the Manna uni- 
versality class (Paczuski and Boettcher, 1996; Chianca et al, 2009). After the defi- 
nition of the model, the discussion below focuses on the model's role in earthquake 
modelling and the attention it received for the spatio-temporal patterns it develops. 

The OFC Model is at home on a two-dimensional square lattice. As in the mod- 
els above, each site ; has a local degree of freedom z,- e R (called the local "pulling 
force"), which is, in contrast to the models above, however, real-valued. As in the 
BTW Model, there are two clearly distinct stages of external driving and internal 
relaxation. During the driving all sites in the system receive the same amount of 
force (sometimes referred to as "continuous" or better "uniform" drive) until one 
site exceeds the threshold i ' = 1, which triggers a relaxation during which no fur- 
ther external driving is applied. In a relaxation or toppling, a site re-distributes a 
fraction of all pulling force evenly among its nearest neighbours which may in turn 
exceed the threshold. The force z, at the toppling site ; is set to and the amount 
arriving at each neighbour is azu where a is the level of conservation. At coordi- 
nation number q, a level conservation less than l/q means that the bulk dynamics is 
dissipative. Boundary sites lose force azi (at corners multiples thereof) to the out- 
side. Because the force re-distributed depends on the amount of pulling force present 
at the site at the time of the re-distribution, the order of updates matters greately, i.e. 
the OFC Model is not Abelian. If a < l/q periodic boundary conditions can be ap- 
plied without losing the possibility of a stationary state, yet normally the boundaries 
are open. The OFC Model is normally initialised by assigning random, independent 
forces from a uniform distribution. 

Sites to topple are identified at the beginning of a timestep and only those have 
been relaxed by the end of it (parallel updates). Unless more than one site exceeds 
the threshold (degenerate maximum) at the beginning of an avalanche, toppling sites 
therefore reside on either of the two next nearest neighbour sublattices of a square 
lattice. 

Again, a separation of time scales is applied, where the relaxation becomes in- 
finitely fast compared to an infinitesimally slow drive. In an actual implementation, 
however, the driving is applied instantaneously and the relaxation takes up most 
(computational time): The driving can be completed in a single step by keeping 
track of the site, say i* with the largest pulling force acting on it. The amount 
of force added throughout the system is thus simply z c — Z;* , triggering the next 
avalanche. 
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Because sweeping the lattice in search of the maximum is computationally very 
costly, 8 the main computational task in the OFC Model is to keep track of the site ex- 
posed to the maximum pulling force. This is a classic computational problem (Cor- 
men et al, 1996), which also is occurs in other models, such as the Bak-Sneppen 
Model (Bak and Sneppen, 1993). The traditional solution is to organise data in a 
tree-like structure and devise methods that allow fast updating and determination of 
the maximum. However, in the OFC Model updating as site's force is much more 
frequent than determination of the maximum and thus a fast algorithm focuses on 
the optimisation of the former at the expense of the latter, i.e. a slightly slower pro- 
cedure to determine the maximum. 

Grassberger (1994) pointed out a number of improvements to a naive, direct im- 
plementation of the OFC Model. Firstly, instead of driving the system uniformly, 
thereby having to sweep the lattice to increase the force on every site by z c — z,* , 
the threshold z c is to be lowered; the amount of force re-distributed at toppling is 
obviously to be adjusted according to the new offset. The second major improve- 
ment Grassberger (1994) suggested was the organisation of forces in "boxes" (some- 
times referred as Grassberger's box-technique), which splits the range of forces 
present in the system in small enough intervals that the search for the maximum 
force succeeds very quickly, yet keeps the computational effort to a minimum when 
re-assigning a box after an update. Other improvements suggested was maintaining 
a stack (Sec. 1.3.1) of active sites, and the use of a scheme to determine neighbour- 
ing sites suitable to the programming language at hand. 

The adjustment of f outlined above has some rather unexpected effects de- 
pending on the numerical precision (Sec. 1.3.3) used in the simulation (Pruessner, 
2012c). As pointed out by Drossel (2002), the OFC Model is extremely sensitive to 
a change of precision; a lower precision seems to enhance or favour phase-locking, 
discussed in the following. 

Most of the studies of SOC models focuses on large-scale statistical features, 
large both in time and space. The analysis of the OFC Model by Socolar et al. 
(1993) Middleton and Tang (1995) and Grassberger (1995) therefore ventured into 
uncharted territory as they studied the evolution towards stationarity in the OFC 
Model on a microscopic scale, analysing the patchy structure of the forces on the 
lattice. 

Firstly, periodic boundary conditions inevitably lead to periodic behaviour in 
time. Below a « 0.18 in a two-dimensional square lattice, (almost) every avalanche 
has size unity. In that extreme case, the period is strictly 1 — qa, because discount- 
ing the external drive, this is the amount of force lost from every site after every site 
has toppled exactly once, as the system goes through one full period. 

The periodicity is broken once open boundaries are introduced. Sites at the edge 
of the lattice have fewer neighbours that charge them, so if every site in the system 
topples precisely once, the force acting on a boundary site is expected to be lower. 
While open boundaries indeed break temporal periodicity, they form, at the same 

8 Not only is the very searching across all sites costly, most of the memory occupied by the lattice 
will not reside in a cache line (as for example most "local" data) and thus has to be fetched through 
a comparatively slow bus. 
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time, seeds for (partially) synchronised patches, which seem to grow from the out- 
side towards the inside, increasing in size. Middleton and Tang (1995) introduced 
the term marginal (phase) locking to describe this phenomenon. 

The temporal periodicity might similarly be broken by introducing inhomo- 
geneities or disorder, effective even at very low levels (Grassberger, 1994; Ceva, 
1995, 1998; Torvund and Fr0yland, 1995; Middleton and Tang, 1995; Mousseau, 
1996). That a spatial inhomogeneity helps forming synchronised patches in space 
can also be attributed to marginal phase locking. 

Because the OFC Model is so sensitive to even the smallest amount of disorder 
and inhomogeneity, its statistics is often taken from very big samples with extremely 
long transients. Many authors also average over the initial state. Drossel (2002) 
suggested that despite these precautions, some of the statistical behaviour allegedly 
displayed by the OFC Model might rather be caused by numerical "noise", also 
a form of inhomogeneity or disorder entering into a simulation. In practise, it is 
difficult to discriminate genuine OFC behaviour from numerical shortcomings and 
one may wonder whether some of these shortcomings are not also present in the 
natural phenomenon the OFC Model is based on. 

That SOC may be applicable in seismology had been suggested by Bak et al. 
(1989b, also Bak and Tang, 1989; Sornette and Sornette, 1989; Ito and Matsuzaki, 
1990) at a very early stage. The breakthrough came with the OFC Model, which is 
based on the Burridge-Knopoff Model of earthquakes (or rather fracturing rocks). 
The latter is more difficult to handle numerically, with a "proper" equation of mo- 
tion taking care of the loading due to spring-like interaction much more carefully. 
The OFC Model, on the other hand, is much easier to update, almost like a cel- 
lular automaton. 9 The context of SOC provided an explanatory framework of the 
scale-free occurrence of earthquakes as described by the Gutenberg-Richter law 
(Gutenberg and Richter, 1954; Olami et al, 1992). Even though exponents both 
in the real-world as well as in the OFC Model seem to lack universality, certain 
scaling concepts, motivated by studies in SOC, have been applied successfully to 
earthquake catalogues (Bak et al, 2002). 

It is fair to say that the OFC Model, to this day, is widely disputed as a bona 
fide model of earthquakes. Its introduction has divided the seismology community, 
possibly because of the apparent disregard of their achievements by the proponents 
of SOC (Bak and Tang, 1989). One of the central claims made initially is that earth- 
quakes are unpredictable if they are "caused" by SOC, which questions the very 
merit of seismology. However, given that SOC is a framework for the understanding 
of natural phenomena on a long time and length scale, providing a mechanism for 
the existence of long temporal correlations, SOC indicates precisely the opposite of 
unpredictability. This point is discussed controversially in the literature to this day 
(Corral, 2003, 2004c,b; Davidsen and Paczuski, 2005; Lindman et al, 2005; Corral 
and Christensen, 2006; Lindman et al, 2006; Werner and Sornette, 2007; Davidsen 
and Paczuski, 2007; Sornette and Werner, 2009). Older reviews (Turcotte, 1993; 
Carlson et al, 1994) help to understand the historical development of the dispute. 



9 Strictly, the OFC Model generally is not a cellular automaton, because the local states n are 
continuous. 
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Hergarten (2002) and more recently Sornette and Werner (2009) have put some of 
the issues in perspective. 

There is not a single set of exponents for the OFC Model, as they are generally 
expected to vary with the level of conservation (Christensen and Olami, 1992). Be- 
cause authors generally do not agree on the precise value of a to focus on, results 
are not easily comparable across studies. Even in the conservative limit, a = l/q, 

little data is available, suggesting z = 1.22(5) 1.253 and D = 3.3(1) 3.01 

(Christensen and Olami, 1992; Christensen and Moloney, 2005). 



1.2 Scaling and numerics 

As a rule of thumb, SOC models are SDIDT systems Jensen (1998): Slowly Driven 
Interaction Dominated Threshold systems. The driving implements a separation of 
time scales and thresholds lead to highly non-linear interaction, which results in 
avalanche-like dynamics, the statistics of which displays scaling, a continuous sym- 
metry. Ideally, the scaling behaviour of an SOC model can be related to some under- 
lying continuous phase transition, which is triggered by the system self-organising 
to the critical point. 

The critical behaviour can be characterised by (supposedly) universal critical 
exponents, the determination of which is the central theme of the present section. 
At the time of the conception of SOC, critical exponents were extracted directly 
from probability density function, (PDFs), often by fitting the data to a straight line 
in double-logarithmic plot. Frequently, such scaling is referred to as "power law 
behaviour". Very much to the detriment of the entire field, some authors restrict 
their research to the question whether an observable displays the desired behaviour, 
without attempting to determine its origin and without considering the consequences 
of such behaviour. Power law behaviour therefore has become, in some areas, a mere 
curiosity. 



1.2.1 Simple scaling 

While studying power laws in PDFs can be instructive, there are far superior meth- 
ods to quantify scaling behaviour. In recent years, most authors have focused on an 
analysis of the moments of the PDFs, as traditionally done in the study of equilib- 
rium statistical mechanics. Not only is this approach more efficient, it also is more 
accurate and mathematically better controlled. Moreover, it is concerned directly 
with an observable (or rather, arithmetic means of its powers), rather than its accu- 
mulated histogram. 

Nevertheless the starting point of a scaling analysis in SOC, is the simple (finite 
size) scaling assumption, 
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9>(s) = as- T &(s/s c ) for s » s , 



(1-3) 



where 3^{s) is the (normalised) probability density function of an observable, s in 
this case, a is a (non-universal) metric factor present to restore dimensional con- 
sistency and accounting for the (microscopic) details of the model, T is a universal 
scaling (or critical) exponent, ^ is a universal scaling function, s c is the upper 
cutoff and sq the lower cutoff. If s is the avalanche size, then T is known as the 
avalanche size exponent, when s is the duration, then x is traditionally replaced by 
a and called the avalanche duration exponent. 

The two cutoffs signal the onset of new physics: Below sq some microscopic 
physics prevails, often a lattice spacing or some other minimal length below which 
discretisation effects take over. Above s c some large finite length scale becomes vis- 
ible, which in SOC is normally controlled by the size of the lattice, so that Eq. (1.3) 
is referred to as finite size scaling. In traditional critical phenomena, s c is controlled 
by the correlation length, beyond which distant parts of the system can be thought 
of as being independent, suggesting the validity of the central limit theorem. 

Strictly, SOC models should always tune themselves to a critical point, so that 
the algebraic, critical behaviour is cut off only by the system size. All scaling in 
SOC therefore is finite size scaling. There are a handful of established SOC mod- 
els, which violate that strict rule, however, such as the Drossel-Schwabl Forest Fire 
Model Drossel and Schwabl (1992a), where an additional parameter has to be tuned 
simultaneously with the system size. 

The physical origin of the scales contained in the metric factor a and the lower 
cutoff so often is the same, yet even with these length scales present, &(s) has an 
arbitrarily wide region where it displays a power-law dependence on s and whose 
width is controlled by s c ; if sq « s « s c , then £P(s) = as~ T+a s~ a %, provided 



Typically, however, a = so that the intermediate region of displays a power 
law dependence with exponent T, which can in principle be extracted as the neg- 
ative slope of 3?{£) in a double logarithmic plot. However, because it is a priori 
unclear whether the scaling function ^(s/s c ) can be approximated sufficiently well 
by a constant "measuring" the exponent T by fitting the intermediate region of a 
double logarithmic plot to a straight line (sometimes referred to as the apparent ex- 
ponent) is very unreliable. If the scaling function displays a power law dependence 
on the argument, a ^ 0, the effective exponent in the intermediate region is T — a. 
One can show that a is non- negative, a > 0, and T = 1 if a > (Christensen et al., 



Figure 1 .3 shows a typical double-logarithmic plot of the PDF in an SOC model. 
The power law region is marked by two dashed lines. The lower cutoff is at around 
io = 50 and the features below that value are expected to be essentially reproduced 
by that model irrespective of its upper cutoff. The spiky structure visible in that 
region is not noise and may, to some extent, be accessible analytically, similar to the 
lattice animals known in percolation (Stauffer and Aharony, 1994). The power law 



(1.4) 
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Fig. 1.3 Example of a double logarithmic plot of the PDF of the avalanche size in an SOC model 
(Data from Pruessner, 2012c). 

region between the two dashed lines can be widened arbitrarily far by increasing 
the upper cutoff s c . Running the same model with increasing s c will reproduce this 
almost straight region beyond which the bump in the data indicates the onset of the 
upper cutoff. 

The upper cutoff in SOC models supposedly depends only on the system size 
and does so in a power-law fashion itself, 

s c {L)=bL D (1.5) 

where b is another metric factor and D is the avalanche dimension. The exponent 
describing the same behaviour for the upper cutoff of the avalanche duration is the 
dynamical exponent z. The four exponents %, D, a and z are those most frequently 
quoted as the result of a numerical study of an SOC model. 

The simple scaling ansatz Eq. (1.3) as well the scaling of the upper cutoff, 
Eq. (1.5), both describe asymptotic behaviour in large s c and L respectively. When 
determining exponents in computer simulations of SOC models, corrections have to 
be taken into account in a systematic manner. While subleading terms are difficult 
to add to the simple scaling ansatz Eq. (1.3), this is routinely done in the case of the 
upper cutoff, 

s c {L) = bL D (l+c l L-'° l +c 2 L- 6>2 ...) (1.6) 

Corrections of this form are referred to as corrections to scaling (Wegner, 1972) 
or confluent singularities. These corrections are discussed further in the context of 
moment analysis, Sec. 1.2.2. 

Although some very successful methods of analysis exist (Clauset et ai, 2009), 
Eq. (1.3) does not lend itself naturally to a systematic quantitative analysis for fixed 
s c . Often, a data collapse is performed in order to demonstrate the consistency of the 
data with simple scaling. According to Eq. (1.3) the PDF 3?{s) for different cutoffs 
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Fig. 1.4 Data collapse of three different data sets similar to the data shown in Figure 1.3. The 
upper cutoff s c is solely controlled by the system size L (Data from Pruessner, 2012c). 



s c produces the same graph by suitable rescaling, in particular by plotting s T 
against x = s/s c , which gives (x). Deviations are expected for small values of s/s c , 
namely for s around so, where Eq. ( 1 . 3) does not apply. Figure 1 .4 shows an example 
of such a collapse using the same model as in Figure 1.3. 

Provided ]im x -^o&(x) = % ^ 0, the region where 3P{s) displays (almost) a 
power law translates into a horizontal, (nearly) constant section in the rescaled plot. 
The graph terminates in a characteristic bump, where the probability density of 
some larger event sizes exceeds that of some large, but smaller ones. This counter- 
intuitive feature is normally interpreted as being caused by system spanning events 
which were terminated prematurely by the boundaries of the system. Had the sys- 
tem been larger, the events would have developed further. In the PDF of a larger 
system thus make up the regular, straight power law region, where the smaller sys- 
tem's PDF displays a bump. Even when the total probability contained in the bump 
is finite but very small, it is enough to account for all events contained beyond it in 
the power law region of an infinite system. 

A data collapse is not unique, as plotting 3^{s) s T f(s/s c ) produces @(x)f(x) 
for any function f(x). In the literature, f(x) is often chosen as f{x) = x~ T so that 
&{s) s r f(s/s c ) = 0P{s)s1. Plotting that data has the fundamental disadvantage that 
s x c usually spans many orders of magnitude more across the ordinate compared 
to 3P{s) s 1 , so that details in the terminal bump are less well resolved. 



1.2.1.1 Binning 



A clean, clear dataset like the one shown in Figure 1.3 is the result of binning. For 
numerical studies of SOC models this is a necessary procedure in order to smoothen 
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otherwise rather rugged histograms. The reason for that ruggedness is the strong 
dependence of the probability density on the event size, with very few large events 
occurring. Because of the power law relationship between event size and frequency, 
their total numbers decrease even on a logarithmic scale. As a result, statistical noise 
visibly takes over, often clearly before the onset of the terminal bump. Statistics for 
large event size is sparse and often little more than a muddle of seemingly unrelated 
data points is visible in the raw data for large events. 

The noise can be reduced by averaging the data for increasingly large event sizes 
over increasingly large "bins", hence the name binning. This is normally done in 
post-processing of the raw data produced in a numerical simulation, by summing 
over all events within a bin and dividing by its size. In principle, the bin sizes could 
be chosen to fit the data; if the bin ranges are [£>,-, £>, + i), then a pure power law 
= as~ x would deposit 



events in each bin i. This number can be made constant by choosing bi = (Bq — 
Biz) 1/(1 ~ T) . Similarly, one might chose the bin boundaries bj "on the fly", i.e. suc- 
cessively increase the bin size until roughly a given number of entries have been 
collected. While those two choices lead to uniformly low statistical errors (assum- 
ing constant correlations), they both suffer from significant shortcomings. Firstly, 
the exponent T to be estimated from the data should not enter into the very prepa- 
ration of the data that is meant to produce the estimate. This problem is mitigated 
by the fact that T may be determined through a separate, independent procedure. 
Secondly and more importantly, both procedures will lead to an increasingly wide 
spacing of data points, which becomes unacceptable towards large event sizes, be- 
cause the abscissa will no longer be defined well enough — if bj + 1 and bi are orders 
of magnitude apart, which s does Eq. (1.7) estimate. Last but not least, to make PDFs 
of different system sizes comparable, the same bi should be used for all datasets. 

The widely accepted method of choice is exponential binning (sometimes also 
referred to as logarithmic binning), where bi = Bo ex P(j3/). Such bins are equally 
spaced on the abscissa of a double logarithmic plot. Because the width of exponen- 
tial bins is proportional 10 to their limits, Eq. (1.7), sparse data can cause a surprising 
artefact, whereby single events spuriously produce a probability density which de- 
cays inversely with the event size, £P(s) ccs~ l , suggesting an exponent of T = 1. A 
typical problem with exponential bins occurs at the small end of the range when 
used for integer valued event sizes, because in that case the bi+\ — bi should not be 
less than 1 . It is then difficult to control the number of bins and thus the resolution 
effectively, because decreasing /3 increases the number of minimally sized bins and 
has highly non-linear knock-on effects on all bin boundaries. The problem is obvi- 
ously much less relevant for non-integer event sizes, such as the avalanche duration. 
However, it is rather confusing to use non-integer bin boundaries for integer valued 
event sizes, because bins may remain empty and the effective bin size cannot be 




(1.7) 



For integer valued bin boundaries, strictly, this holds only approximately. 
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derived from bj + \ — b,. For example a bin spanning b( + \ — b, = 0.9 may not contain 
a single integer, whereas bj + \ — bj = 1.1 may contain two. 

It is obviously advantageous to perform as much as possible of the data manipu- 
lation as post-processing of raw simulation data. Efficiency and memory limitations, 
however, normally require a certain level of binning at the simulation stage. When 
event sizes and frequencies spread over 10 orders of magnitude a simple line of 
code 11 

histogram[size]++; /* one count for size in the histogram */ 

would require histogram to have a precision of more than 32 bits. Normally such 
counters are implemented as integers, which would need to be a long long intin 
the present case. The memory required for 10 10 of these 64 bit numbers (about 75 
GB) exceeds by far the memory typically available in computers in common use at 
the time of writing this text (2012). Writing every event size in a list, eventually to 
be stored in a file, is rarely an alternative, again because of the enormous memory 
requirements and because of the significant amount of computational time post- 
processing would take. 

Consequently, some form of binning must take place at the time of the simula- 
tion. In principle, any sophisticated binning method as used during post-processing 
can be deployed within the simulation, yet the risk of coding errors ruining the 
final result and the computational effort renders this approach unfeasible. The es- 
tablished view that complicated floating point operations such as log or pow are too 
expensive to be used regularly in the course of a numerical simulation has experi- 
enced some revision over the last decade or so, as techniques like hyperthreading 
and out-of-order execution are commonly used even in the FPU. Nevertheless, inte- 
ger manipulation, often doable within a single CPU cycle, remains computationally 
superior compared to floating point manipulation. Even some of the rather archaic 
rules remain valid, such as multiplications being computationally more efficient than 
divisions, as they can be performed within a short, fixed number of cycles. Further 
details can be found in the appendix at the end of the chapter. 



1.2.2 Moment analysis 

By far the most powerful technique to extract universal features of an SOC model 
is a moment analysis (De Menech et al, 1998). Traditionally, the numerical in- 
vestigation of critical phenomena has focused on moments much more than on the 
underlying PDF, even when the former are often seen as the "result" of the latter. 
Mathematically, no such primacy exists and one can be derived from the other un- 
der rather general conditions (Feller, 1966, Carleman's theorem in). In general one 
expects that a finite system produces only finite event sizes, i.e. that finite systems 
have a sharp cutoff of the "largest possible event size". While very physical, this rule 

11 All explicit examples in this chapter are written in C, which is the most widely used program- 
ming language for numerical simulations, as long as they are not based on historic Fortran code. 
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finds its exception in residence times, when particles get "buried" in a "pile" over 
long periods. In the Oslo Model, some of these waiting time distributions seem to 
be moderated by scaling functions that are themselves power laws and may possess 
upper cutoffs exponential in the system size (Dhar and Pradhan, 2004; Pradhan and 
Dhar, 2006, 2007; Dhar, 2006). 

Assuming, however, that all moments 



J '00 
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(1.8) 

Jo 

exist, i.e. are finite, then for n + 1 > x 

(s")^as l c + "^g„ (1.9) 

where ~ is used to indicate equivalence to leading order in large s c . Moments with 
n < T — 1 are not dominated by the scaling in s c , i.e. they are convergent in large s c . 
The (asymptotic) amplitudes g n are defined as 

g n = x"~ T Sf(x) (1.10) 
Jo 

expected to be finite for all n > 0. There is an unfortunate confusion in the literature 
about the (spurious) consequences of (s°y = 1 scaling like s\~ x gQ. If T > 1, then 
the leading order of <^s°) is not given by Eq. (1.9). 

The only scaling in SOC is finite size scaling, i.e. the upper cutoff is expected to 
diverge with the system size, Eq. (1.5), so that moments scale like 

(s")~ab l+ '- T L D{1+n -^g n . (1.11) 

Neither a nor b are universal and neither are the g n unless one fixes some features of 
(x) such as its normalisation and its maximum. To extract universal characteristics 
of (x), moment ratios can be taken for example 

gn _ lgn+1 

+ corrections (1-12) 



or 



(s n )(s)"- 2 ~ n ~ 2 



g n + corrections , (1-13) 



which is particularly convenient because of its very simple form when fixing g\ = 
g2 = 1 by choosing the metric factors a and b appropriately. 

The most important result of a moment analysis, however, are the universal expo- 
nents D and T and corresponding pairs for avalanche duration (z and a respectively), 
as well as the area (normally D a and T a ) etc.. This is done in a three step process. 
Firstly, the SOC model is run with different systems sizes, typically spaced by a 
factor 2, or 2,5,10. It can pay to use slightly "incommensurate" system sizes to 
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identify systematic effects, for example due to boundary effects being particularly 
pronounced in system sizes that are powers of 2. A typical simulation campaign 
would encompass 10 to 15 system sizes, of which maybe only 6 to 10, stretching 
over two to four orders of magnitude 12 will be used to produce estimates of expo- 
nents. The result of the simulation are estimates for the moments of the relevant 
observables together with their error (see below). 

Secondly, the moments of the event sizes distribution, (s n \ are fitted against a 
power law in L (which is the parameter controlling s c ) with corrections, 

{s n )=A L^+A i L 11 "- 031 + ... (1.14) 

with positive exponents known as confluent singularities; in particular jx n — a>\ is 
sometimes referred to as a sub-dominant exponent. The introduction of such correc- 
tions to scaling goes back to Wegner (1972), who applied it in equilibrium critical 
phenomena. The Levenberg-Marquardt algorithm (Press et al, 2007) is probably 
the fitting routine most frequently employed for matching the estimates (with their 
error bars) from the simulation to the fitting function Eq. (1.14). There are a number 
of problems that can occur: 

• Unless the result is purely qualitative, a good quality fit cannot be achieved with- 
out good quality numerical data, that includes a solid estimate of the numerical 
error, i.e. the estimated standard deviation of the estimate. 

• The very setup of fitting function Eq. (1.14) (sometimes referred to as "the 
model") can introduce a systematic error; after all it is only a hypothesis. 

• If n > x — 1 is very small, corrections due to the presence of the lower cutoff (so, 
Eq. (1.3)) can be very pronounced. 

• The error stated for the fitted exponents alone can be misleading. If Eq. (1.14) is 
very constraining, the error will be low, but so will the goodness-of-fit. 

• Too many fitting parameters allow for a very good goodness of fit, but also pro- 
duce very large estimated statistical errors for the exponents. 

• Fitting against a function with many parameters often is highly dependent on the 
initial guess. In order to achieve good convergence and systematic, controlled 
results, it may pay off to fit the data against Eq. (1.14) step-by-step, using the 
estimates obtained in a fit with fewer corrections as initial guesses for a fit with 
more corrections. 

• In most cases, there is little point in having as many parameters as there are 
data points, as it often produces a seemingly perfect fit (goodness-of-fit of unity), 
independent of the input data. 

• Extremely accurate data, i.e. estimates for the moments with very small error 
bars, may require a large number of correction terms. 



One might challenge the rule of thumb of the linear system size L having to span at least three 
orders of magnitude — in higher dimensions, say d = 5, spanning three orders of magnitude in lin- 
ear extent leads to 15 orders of magnitude in volume, which might be the more suitable parameter 
to fit against. 
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• It can be difficult to force the corrections a), to be positive. It is not uncommon to 
fix them at certain reasonable values such as 0), = i or a); = i/2. Alternatively, they 
can be introduced differently, writing them, for example, in the form a); = i + \ 6h\ . 

• If finite size scaling applies, the relative statistical error for any moment scales 
like (s 2 ") / (s") 2 <xL D ( T ~ l \ assuming that a 2 (s n ) scales like (^ 2 "), which it cer- 
tainly does for z > 1. At z = 1 the scaling of o 2 (s n ) may be slower than that of 
(s") 2 . While L D( - T -V does not depend on n, the amplitude of the moments does, 
leading normally to an increase of the relative error with n. 

In some models the first moment of the avalanche size displays anti-correlations 
and thus faster numerical convergence as found in a mutually independent sample 
(Welinder et al, 2007). In many models, the average avalanche size (s) is known 
exactly, in one dimension often including the confluent singularities (Pruessner, 
2012a). These exact results can provide a test for convergence in numerics and also 
provide a scaling relation 

D(2-t) = jUi (1.15) 

If ji\ is known exactly = 2 for bulk driving Manna, Oslo and Abelian BTW 
Models, fi\ = 1 for boundary drive), then Eq. (1.15) gives rise to a scaling relation. 
Normally, there are no further, strict scaling relations. However, the assumption of 
narrow joint distributions suggests D(z — 1) = z(ce — 1) etc. (Christensen et al, 
1991; Chessa et al., 1999a; Pruessner and Jensen, 2004). If the exponent [i\ is given 
by a mathematical identity and (s) serves as an analytically known reference in the 
numerical simulation, then jii should not feature in the numerical analysis to extract 
the scaling exponents D and z. Rather, when fitting n„ versus D(l + n — z), this 
should be replaced by D(n — 1 ) + jii . 

Fitting ji„ in a linear fit (without corrections) against D(l + n — z) (or against 
D{n — 1) + jUi if jii is known exactly) is, in fact, the third step in the procedure 
described in this section. In principle, the n > % — 1 do not need to be integer valued. 
They have to be large enough to avoid a significant corrections due to the lower 
cutoff, and small enough to keep the relative statistical error small. Non-integer « 
can be computationally expensive, as they normally require at least one library call 

of pow. 

While each estimate jx n for every n should be based on the entire ensemble, 
considering them together in the same fit to extract the exponents D and z introduces 
correlations, which are very often unaccounted for. As a result both goodness-of-fit 
as well as the statistical error for the exponents extracted are (unrealistically) small. 

There are a number of strategies to address this problem. The simplest is to up- 
scale the error of the fi„ as if every estimate was based on a separate, independent 
set of raw data. Considering M moments simultaneously, their error therefore has to 
scaled up by a factor \[M (Huynh et al., 2011). In a more sophisticated approach, 
one may extract estimates from a series of sub-samples (Efron, 1982; Berg, 1992, 
2004). 

It often pays to go through the process of extracting the exponents D and z at 
an early stage of a simulation campaign, to identify potential problems in the data. 
Typical problems to watch out for include 
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• Corrections are too strong because system sizes are too small. 

• Results are too noisy because sample sizes are too small, often because the sys- 
tem sizes are too big. 

• Results have so little noise that fitting functions need to contain too many free 
parameters. 

• Too few data points (i.e. too few different system sizes L or different moments 
n). 

• Large event sizes suffer from integer overflow, resulting in seemingly negative or 
very small event sizes. 

• Data identical in supposedly different runs, because of using the same seed for 
the random number generator. 

• Transients chosen too short. 



1.2.3 Statistical errors from chunks 

One of the key-ingredients in the procedures described above is a reliable estimate 
for the statistical error of the estimates of the individual moments. The traditional 
approach is to estimate the variance, a 2 (s") = (s 2 "^ — (s") 2 of each moment, so 
that the statistical error of the estimate of (s") is estimated by a 2 \s n ) / 'a/ 'N / '(2t + 1), 
where Nj (2t + 1) is the number of effectively independent elements in the sample 
with correlation time t. 

This approach has a number of significant drawbacks. Firstly, each moment <V') 
requires a second moment, <(s 2n ), to be estimated as well. Considering a range of 
moments, this might (almost) double the computational effort. Rather dissatisfy- 
ingly, the highest moment estimated itself cannot be used to extract its finite size 
scaling exponent fi„, because its variance is not estimated. Furthermore, because 
of their very high powers, the moments entering the estimates of the variances and 
thus the variances themselves have large statistical errors and are prone to integer 
overflow. 

Estimating the effective number of independent elements in the sample is a hur- 
dle that can be very difficult to overcome. Usually, it is based on an estimate of the 
con-elation time t. If (s/sy) — (s) 2 = c 2 (s)exp(— |;' — j\/z), then the variance of the 



estimator 




(1.16) 



of (s) for N » T is 



(j2 ^ = pE(<^>-< i > 2 )* 



JV(l-exp(-l/T)) 



1 +exp(-l/r) 




2t+1 



a 2 (s) (1.17) 



as if the sample contained only N/ (2t + 1 ) independent elements. 
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The main difficulty of this strategy is a reliable estimate of % which often cannot 
be easily extracted from (siSj) — (s) 2 because of noise and the presence of other 
exponential contributions, of which exp(— \i — is the slowest decaying one. 
Moreover, in principle T has to be measured for each observable separately (even 
when it makes physically most sense to assume that the system is characterised by 
a single correlation time). 

To avoid these difficulties, one may resort to a simple sub-sampling plan. As 
discussed below (also Sec. 1.3.5), it is a matter of mere convenience and efficiency 
to repeatedly write estimates of moments based on a comparatively small sample 
into the output stream of a simulation and reset the cumulating variables. In the 
following these raw estimates based on a small sample are referred to as chunks. If 
their sample size is significantly larger than the correlation time, then each of these 
estimates can be considered as independent and the overall estimates based on it has 
its statistical error estimated accordingly. For example, if m, with ; = 1,2, ... ,M are 
estimates of <V') all based on samples of the same size N, say m; = YTj ' s 1j w ^ m S U 



the jth element of the i sample, then the overall unbiased and consistent estimator 
(Brandt, 1998) of <» is 



One crucial assumption above is that the m, are independent, which can always be 
achieved by merging samples. As long as M remains sufficiently large, one may be 
generous with the (effective) size of the individual samples (Flyvbjerg and Petersen, 
1989). 

Chunks also allow a more flexible approach to determining and discarding tran- 
sient behaviour from the sample supposedly taken in the stationary state. The tran- 
sient can be determined as a (generous) multiple of the time after which (ideally all 
or several) observables no longer deviate more from the asymptotic or long time av- 
erage than their characteristic variance. Where observables are known exactly (e.g. 
the average avalanche size Pruessner, 2012a), they can be used as a suitable refer- 
ence. Figure 1.5 shows the transient behaviour of the average avalanche size in a 
realisation of the Manna Model. A more cautious strategy is to consider a series of 
different transients and study the change in the final estimates (with their estimated 
error) as a function of the transient discarded. 




(1.18) 



which has an estimated standard deviation of (m 2 



m 2 )/(M— 1) where 




(1-19) 



32 



messner 




Fig. 1.5 Example of the transient behaviour of an observable (here the average avalanche size in 
the one-dimensional Manna Model with L = 65536) as a function of the chunk index in a log-lin 
plot (data from Huynh et al., 201 1). The straight dashed line shows the exact expected average (s), 
Eq. (1.1). The arrow indicates the chunk from where on stationarity is roughly reached. A generous 
multiple of that time should be taken as the number of chunks to discard in order to ensure that 
correlations (and thus dependence on the initial setup) are essentially overcome. 

1.3 Algorithms and data organisation 

In the following, a range of numerical and computational procedures are discussed 
that are commonly used in the numerical implementation of SOC models (for a 
more extensive review see Pruessner, 2012c). Some of them are a matter of common 
sense and should be part of the coding repertoire of every computational physicist. 
However, it is not always entirely obvious how these "standard tricks" are used for 
SOC models. 

In the following, the focus is on computational performance, which often comes 
with the price of lower maintainability of the code. The amount of real time spent 
on writing code and gained by making it efficient, should account for the time spent 
on debugging and maintaining it. 

Most of the discussion below is limited to algorithmic improvements. The aim 
is produce code that communicates only minimally with the "outside world", be- 
cause in general, interaction with the operating system, as required for writing to a 
file, is computationally expensive and extremely slow. The UN*X operating system 
family (including, say, Linux and Mac OS X) distinguishes two different "modes" 
by which an executable keeps the CPU busy: By spending time on the (operating) 
system and by spending it in "user mode". Roughly speaking, the former accounts 
for any interaction between processes, with external controls or peripherals, includ- 
ing writing files. The latter accounts for the computation that takes place solely on 
the CPU (ALU, FPU, GPU, etc.) and the attached RAM. Tools like time and li- 
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brary functions like getrusage provide an interface to assess the amount of various 
resources used, while being themselves or resulting in systems calls. 

Of course, the literature of computational physics in general is vast. Reviews and 
texts that are of particular use in the present context include Kernighan and Ritchie 
(1988); Cormen et al. (1996); Knuth (1997); Newman and Barkema (1999); Berg 
(2004); Landau and Binder (2005); Press et al. (2007). 



1.3.1 Stacks 

The definition of most SOC models makes no reference to the method to identify 
active sites, i.e. sites that are due to be updated. In principle, an implementation of 
an SOC model could therefore repeatedly scan the entire lattice to find the relevant 
sites. This is, however, very inefficient and therefore should be avoided. Instead, the 
index of active sites (or their coordinates) should be organised in a list. Every site 
in that list is subsequently updated. Moreover, it is often very important to know 
whether a site is maintained in the list or not. Sometimes this can be determined 
implicitly (for example, when a site is guaranteed to reside on the list from the mo- 
ment its height exceeds the threshold), sometimes this is done explicitly by means 
of a flag associated with the site. The following contains a more detailed discussion 
of the various techniques available. 

The most commonly used form of a list is a stack, called so, because this is 
how it appears to be organised. It consists of a vector, say int stack [stack_size 
], of pre-defined size stack_size. It must be large enough to accommodate the 
maximum number of simultaneously active sites. Simulating large lattices, a balance 
has to be struck between what is theoretically possible and what is happening in 
practise. 

The type of the stack, vector of int in the example above, is determined by the 
type it is meant to hold. If it holds the index of active sites, it is likely to be be 
int, but it may also hold more complex objects, say, coordinates of active particles 
(but see below). The number of objects currently held by the stack is stored in int 

stack_height . 

If stack_size is smaller than the theoretical maximum of active sites, int 
stack_height has to be monitored as to prevent it from exceeding stack_size 
. The outcome of the simulation is undefined if that happens, because the ex- 
act position in memory of stack [stack_size] is a priori unknown. If therefore 
stack_height exceeds stack_size, memory has to be extended one way or an- 
other. For example, one may use reaiioc ( ) , which assumes, however, that enough 
memory is actually available. Modern operating systems all provide virtual memory 
which is transparently supplemented by a swap file residing on the (comparatively 
slow) hard drive. This is to be avoided because of the computational costs associ- 
ated. It may thus pay off for the process itself to make use, temporarily, of a file to 
store active sites. The alternative to abandon the particular realisation of the simula- 
tion introduces a bias away from rare events which is likely to have significant effect 
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on observables. The same applies obviously if activity is suppressed if it reaches the 
maximum level. 

There are two fundamental operations defined on a stack, 

#define PUSH(a) stack[stack_height++]=(a) 
#define POP(a) (a)=stack[ stackjieight] 

where push (a) places (a) on the stack and pop takes an element off. The under- 
lying idea is literally that of a stack: When a site becomes active, its index goes on 
a pile (push) so that each index number on that pile represents a site waiting to be 
updated. When that happens, it is removed from the pile (pop). 

It simplifies the code greatly if all objects on the stack are, in a sense, equivalent. 
For example, all sites on a stack are active. Guaranteeing this is not necessarily 
trivial, because the manipulation of one item on the stack may affect the state (and 
thus the eligibility) of another item on the stack. It is therefore advisable to ensure 
that all elements on the stack are distinct. In SOC models that means that active sites 
enter the stack exactly once, namely when the turn active. If an active site is charged 
again by a toppling neighbour, a new copy of its index is not placed on the stack. In 
the Manna Model, for instance, the single line of code to place objects on the stack 
could be 

if (z[i]++==l) {PUSH(i);} 

so that the index i of a site enters when it is charged while its height z is at the 
critical value z c ■ The line should not read if (z [i] ++>=l) push (i) ; . 

Unfortunately, the very data structure of a stack, which in the present context may 
better be called a LIFO (last in, first out), suggests a particular procedure to explore 
active sites, namely a depth first search (DFS); Whenever a toppling site activates 
its neighbours, one of them will be taken off first by the next call of pop, toppling 
in turn. Activity thus spreads very far very quickly, then returning, then spreading 
far again, rather than "burning locally". In fact, in DS-FFM a DFS is probably the 
simplest way of exploring a cluster of trees. 

The alternative, a breadth first search (BFS) requires slightly greater computa- 
tional effort because it normally makes use of a FIFO (first in, first out). The last 
object to arrive on a FIFO is the last one to be taken off, exactly the opposite or- 
der compared to a stack. Naively, this may be implemented by removing items 
from the front, stack [ ] , and using memmove 13 to feed it from the end, lowering 
stack_height. This approach, however, is computationally comparatively costly. 
A faster approach is to organise the stack in a queue, organised in a ring (circular 
buffer) to keep it finite, where a string of valid data grows at the end while retreating 
from the front. 

In Abelian models, where the statistics of static features of avalanches, such as 
size and area, do not depend on the details of the microscopic dynamics 14 , working 



Dedicated library functions like memmove and memcpy are generally much faster than naive 
procedures based on loops, although the latter can be subject to significant optimisation by the 
compiler. 

14 But note the strict definition of Abelianness discussed on p. 7. 
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through the stack using pop may be acceptable. Where temporal features are of 
interest too, the microscopic dynamics must implement a suitable microscopic time 
scale. Often the microscopic timescale is given by Poissonian updates, for example 
by active sites toppling with a Poissonian unit rate. 

In principle that means that waiting times between events (sites toppling) are 
themselves random variables. If a faithful representation of the microscopic time 
is desired, then the random waiting times can be generated by taking the negative 
logarithm of a random number drawn from a uniform distribution on (0, 1]. If an 
approximate representation of the Poisson processes is acceptable (which, in fact 
converges to the exact behaviour in the limit of large numbers of active sites, see 
Liggett, 2005), then elements are taken off the stack at random and time is made to 
progress in steps of l . / stack_height . If stackjieight remains roughly constant, 
than on average stack_height events occur per unit time as expected in a Poisson 
process. A simple implementation reads 

int rs_pos; 

#define RANDOM_POP(a) rs_pos=rand() % stackjieight; (a)=stack[rs_pos]; POP(stack[ 
U rs.pos]) 

where the last operation, pop (stack [rs_pos] ) overwrites the content of stack [ 
rs_pos] by stack [stack_height-i] decrementing stack_height at the same 
time. When selecting the random position on the stack via rs_pos=rand ( ) % 
stack_height a random number generator has to be used (Sec. 1.3.4), which only 
for illustrative purposes is called rand ( ) here. 

One consequence of the constraint of distinct objects on the stack is that a site 
may need to topple several times before being allowed to leave the stack. In Abelian 
models some authors circumvent that by placing a copy of the site index on the stack 
every time a pair of particles has to be toppled from it, which can be implemented 
easily by removing an appropriate number of particles from the site each time it 
enters the stack. As a result, however, stacks may become much larger, i.e. a greater 
amount of memory has to be allocated to accommodate them. 

Depending on the details of the microscopic dynamics, an possible alternative is 
to relax a site completely after it has been taken off the stack, for example in the 
Manna Model: 

while (stackjieight) { 
RANDOM_POP(i); 
do{ 

topple(i);/* Site i topples, removing two particles from i. */ 
avalanche_size++; /* avalanchejiize counts the number of topplings. */ 
} while (z[i]>l); 

} 

where topple (i) reduces z [i] by 2 each time. If the avalanche size counts the 
number of topplings performed, avaianche_size has to be incremented within the 
loop. Counting only complete relaxations would spoil the correspondence with exact 
results. 
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An alternative approach with different microscopic time scale is to topple a site 
on the stack only once, and take it off only once it is fully relaxed. This approach 
requires some "tempering" with the stack: 

while (stackjieight) { 
i=rand() % stackjieight; 
topple(stack[i]); 
if (z[i]<=l) POP(stack[i]); 

} 

In systems with parallel update, where all sites at the beginning of a time step 
have to be updated concurrently before updating the generation of sites that have 
been newly activated, a red-black approach (Dowd and Severance, 1998) can be 
adopted. This requires the use of two stacks, which have to be swapped after com- 
pleting one: 

int *stack, stack_height=0; 

int rb.stack[2][STACK_SIZE], next_stack_height; 
int current.stack, next_stack; 

#define NEXT_PUSH(a) rb_stack[next_stack] [next_stack_height++]=(a) 
#define NEXT_POP(a) (a)=rb_stack[next_stack][ next_stack_height] 

current_stack=0; 
next_stack= 1 ; 

stack=rb_stack[current_stack] ; 

PUSH(i); 

for (;;) { 
while (stackjieight) { 

POP(i); 

NEXT_PUSH(j); 
} "' 

if (next_stackjieight==0) break; 

/* Swap stacks. */ 

stack Jieight=next_stackJieight; 

next_stackJieight=0; 

current_stack=next_stack; 

stack=rb_stack[current_stack] ; 

next_stack= 1 — next.stack; 

} 

/* Both stacks are empty. */ 

The use of the pointer stack is solely for being able to use the macros push 
and pop defined earlier. Otherwise, it might be more suitable to define macros 
Cijrrent_push and current_pop corresponding to next_push and next_pop. 

A stack should also be used when determining the area of an avalanche, i.e. the 
number of distinct sites toppled (or visited, i.e. charged). To mark each site that has 
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toppled during an avalanche and to avoid double counting, a flag has to be set, say 
visited [i] =1 or site [ i ] . visited=i (see Sec. 1.3.2). Counting how often the 
flag has been newly visited then gives the avalanche area. However, in preparation 
for the next avalanche, the flags have to be reset. This is when a stack comes handy, 
say 

int area_stack[SYSTEM_SIZE]; 
int area_stack_height=0; 

#define AREA_PUSH(a) area_stack[area_stack_height++]=(a) 
#define AREA_POP(a) (a)=area_stack[ area_stack_height] 

/* For each toppling site. */ 
if (visited[i]==0) { 

visited[i]=l; 

AREA_PUSH(i); 

} 

/* After the avalanche has terminated. 
* areajitackjieight is the avalanche area. */ 

/* Re— initialise */ 

while (area_stack_height) { 

AREA_POP(i); 

visited[i]=0; 

} 

In the example above, the area is tracked implicitly in area_stack_height . The re- 
initialisation can be further improved using while (area_stack_height ) visited 
[area_stack [ — area_stack_height ] ] =0. 



1.3.2 Sites and Neighbours 

In SOC models, every site has a number of properties, most importantly the local 
degree of freedom, but also (statistical) observables which are being measured and 
updated as the simulation progresses. Other information associated with each site 
are flags (such as the one mentioned above to indicate whether a site had been 
visited) and even the neighbourhood (discussed below). In fact, the site itself may 
be seen as the key associated with all that information. That key might represent 
information in its own right, say, the coordinate, it might be an index of a vector, or 
a pointer. 

1.3.2.1 Pointers and structures 



A word of caution is in order with regard to pointers. The programming language C 
lends itself naturally to the use of pointers. However, code on the basis of pointers 
is difficult to optimise automatically at compile time. Depending on the quality of 
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the compiler and the coding an index based implementation (which is also more 
portable) may thus results in faster code than the seemingly more sophisticated im- 
plementation based on pointers. 

That said, in theory placing pointers on the stack, which gives immediately access 
to a relevant object should be faster than using indices, which are effectively an 
offset relative to abase: b=z [ stack [i] ] might result in machine code of the form b 
=* ( z + * (stack+i) ) which contains one more addition than b=*stack [i] resulting 
in b=** (stack+i) if stack is a vector of pointers. 

Similar considerations enter when using structures, which provide very conve- 
nient and efficient ways of organising and encapsulating data associated with each 
site. For example 

struct site_struct { 
int height; 
char visited; 

}; 

defines a structure with two members, height and visited. Declaring a variable 
struct site_struct site [10] allows the individual elements to be accessed 
in a structured way, say site[i] .height++, site[i] .visited=l. There are a 
number of computational drawbacks, which are, however, normally outweighed by 
the better maintainability of the code. 

• Depending on the platform and the compiler, padding might become necessary, 
i.e. some empty space is added to the structure (Sec. 1.3.2.2, p. 39). The memory 
requirements of the structure is thus greater than the memory requirements for 
each variable when defined individually. 

• Again depending on the platform as well as the compiler, without padding some 
operations on some types may require more CPU cycles (in particular when float- 
ing point types are used). 

• Members within the structure are accessed similar to elements in a vector, namely 
by adding an offset. Access to the first member (where no offset is needed, site 
[i] .height in the example above) can thus be faster than access to the other 
members (site[i] .visited above). Because of that additional addition, the 
approach is often slower than using separate vectors for each member of the 
structure. 



1.3.2.2 Neighbourhood information 

It can be convenient, in particular for complicated topologies or when the neighbour- 
hood information is computed or supplied externally, to store information about the 
local neighbourhood in a site structure, for example: 

struct site_struct { 

int neighbour[MAX_NEIGHBOURS]; 
int numjieighbours; 

}; 
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Because of the significant memory requirements, this is often not viable for large 
lattices. Again, instead of addressing neighbours by their index, pointers can be 
used, which often produces very efficient and elegant code. 

The neighbours of each site thus are calculated and stored at the site only once. 
The strategy of pre-calculated neighbourhoods goes back to the very beginning of 
computational physics, when access to memory was much faster than doing such 
calculations on-the-fly. 15 This, however, has changed. It can be much faster to deter- 
mine a neighbourhood on-the-fly than looking it up, unless, of course, the topology 
is so complicated that it becomes computationally too costly. Unfortunately, it is 
often difficult to try out different implementations (lookup tables and calculation on 
the fly), as the setup of a neighbourhood is at the heart of a lattice simulation. 

As for calculating neighbourhoods, in one dimension the index of a site, which is 
strictly only a key to access the information, is often associated with its position on 
a one-dimensional lattice. Actual computation takes place only at boundaries. If the 
right neighbour of site i in the bulk is i+l, it may not exist on the right boundary or 
be if periodic boundary conditions (PBC) apply in an implementation in C where 
the index of a vector of size length can take values from to length- l. Simi- 
larly, the left neighbour is i-i in the bulk and length-i at i=0 in case of periodic 
boundaries. Those are most easily implemented in the form left= (i+LENGTH-i) % 
length and right= (i+l ) %length respectively using a modulo operation. The shift 
by length in the former avoids problems with negative indices at i=0. 

A less elegant but often faster implementation is to determine whether a site is at 
the boundary before assigning the value for the neighbour, such as 

if (i==0) left=LENGTH— 1 ; 
else left=i— 1; 

or just ieft= (i==0) ?length-i : i-i, which is more readable. This method is 
also more flexible with respect to the boundary condition implemented. Reflecting 
boundary conditions, for example are implemented by ieft=(i==0) ? l : i-l. 
Open boundary conditions, on the other hand, might require special attention. If 
possible, they are best implemented using padding, i.e. by pretending that a neigh- 
bouring site exists, which, however, cannot interact with the rest of the lattice, for 
example, by making sure that it never fulfils the criterion to enter the stack. Such a 
site may need to be "flushed" occasionally to prevent it, for example, from fulfilling 
the criterion due to integer overflow. One might either assign one special site, say 
the variable dump in ieft= (i==0 ) ? dump : i-l or allocate memory for length+2 
sites with an index from to length+i, with valid sites ranging from l to length 
with sites and length +1 receiving charges without toppling in turn. This proce- 
dure also allows a very efficient way to determine the number of particles leaving 
the system, the drop number (Kadanoff et al., 1989). 

Usually only in higher dimensions, one distinguishes reflecting boundary con- 
ditions, where the particle offloaded is moved to another site (normally the mirror 
image of the "missing" site), and "closed" boundary conditions, where the number 



Back in the days when lookup tables for modulo operations were in fashion. 
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of nearest neighbours is reduced and shed particles are evenly re-distributed among 
them. 

Most of the above techniques remain valid in higher dimension, where the data 
can be organised in either a one-dimensional vector or a multidimensional vector. 
The former strategy makes use of macros of the form 

#define COORDINATE2INDEX(x,y,z) ((x)+(LENGTH_X*((y)+LENGTH_Y*(z)))) 
#define INDEX2COORDINATE(i,x,y,z) z=(i)/(LENGTH_X*LENGTH_Y),y=((i)/ 
U LENGTH_X)%LENGTH_Y,x=(i)%LENGTH_X 

The use of the coma operator in the second macro helps to avoid errors when 
omitting curly brackets in expressions like if (l) index2coordinate (i,x,y, z 
) ; . Where stacks are used to hold coordinates, the multiple assignments needed to 
store and fetch all of them may computationally outweigh the benefit of not having 
to calculate coordinates based on a single index. 

The two biggest problem with the use of multi-dimensional vectors is their am- 
biguity when used with fewer indices and the consistency when passing them to 
functions. Both subtleties arise because of the logical difference between a vector of 
pointers to a type and the interpretation of a lower-dimensional variant of a multi- 
dimensional vector. While C makes that distinction, there is no syntactical difference 
between the two. For example 

inta[2][10]; 
a[0][5]=7; 

is a multi-dimensional vector using up 2*i0*sizeof (int) sequential bytes of 
memory. Each a [ i ] is the starting address of each row i = 0, 1 . On the other hand 

int *a[2]; 

introwl[10], row2[10]; 
a[0]=rowl; a[l]=row2; 

a[0][5]=7; 

makes a a vector of pointers, using up 2*sizeof (* int) bytes of memory, while 
each row uses I0*sizeof (int) bytes. Both snippets of code declare a to be com- 
pletely different objects, yet, for all intents and purposes in both cases a will be- 
have like a two-dimensional array. That is, until it is to be passed as an argument 
to another function. In the first case, that function can be declared by function ( 
int array [2] [10] ) , informing it about the dimensions of the array, and subse- 
quently called using function (a) . The two-dimensional vector a will behave as in 
the calling function. In fact, the function will even accept any other vector, lower 
dimensional or not, passed on to it as an argument (even when the compiler may 
complain). 

In the second case, a is a vector of pointers to int, and so a function taking it as 
an argument must be declared in the form function (int **a) , using additional 
arguments or global constants (or variables) to inform it about the size of the vec- 
tor. The two versions of the functions are incompatible, because a two-dimensional 
vector is really a one-dimensional vector with a particularly convenient way of ad- 
dressing its components. In particular, the two-dimensional vector cannot be passed 
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to the function designed for the second case using, say, function (&a) or function 
( (int **) a) . 

While these issues normally are resolved at the time of coding they can cause 
considerable problems when the memory allocation mechanism for the vector is 
changed. This happens, in particular, when lattice sizes are increased during the 
course of a simulation campaign. Initially, one might be tempted to define a lattice 
globally (stored in BSS or data segment) or as automatic variables taken from the 
stack, choosing a multi-dimensional array for convenience. Later on, they make be 
taken from the (usually much bigger) heap using maiioc ( ) , at which point the way 
they are accessed may have to be changed. The latter approach is the most flexible 
but possibly not the most convenient way of allocating memory for large items. 

Finally, it is advisable to scan sites (when sweeping the lattice is unavoidable or 
scanning through a local neighbourhood) in a way that is local in memory and thus 
cache. The first option, declaring a two-dimensional vector in a single step, makes 
that more feasible than the second option, where different rows might end up at very 
different regions of memory. Not using higher dimensional vectors at all, however, 
is probably the best performing option. 



1.3.3 Floating Point Precision 

Very little and at times too little attention is being paid to the effect of limited float- 
ing point precision. Most SOC models can be implemented fully in integers even 
when their degrees of freedom are meant to be real valued, such as the Zhang Model 
(Zhang, 1989), the Bak-Sneppen Model (Bak and Sneppen, 1993) or the Olami- 
Feder-Christensen Model Model (Olami et ah, 1992). In case of the latter, floating 
point precision has been found to significantly affect the results (Drossel, 2002). 

Where random floating point numbers are drawn, they might in fact contain 
much fewer random bits than suggested by the size of their mantissa. In that case, 
an implementation in integers is often not only faster but also "more honest". Where 
rescaling of variables cannot be avoided and occurs frequently, multiplying by a 
constant inverse often produces faster code than division. 

Over the last decade or so, the floating point capabilities of most common CPUs 
have improved so much, however, that the difference in computational costs be- 
tween integers and floating point arithmetics is either negligible or not clear-cut. 
The most significant disadvantage of the latter is the limited control of precision 
that is available on many platforms. 

The levels of precision as defined in the IEEE standard 754 that are very widely 
used are single, double and extended. They refer to the number of bits in the man- 
tissa determined when floating point operations are executed, i.e. they are the preci- 
sion of the floating point unit (FPU). The precision the FPU is running at depends 
on platform, environment, compiler, compiler switches and the program itself. Some 
operating systems offer an IEEE interface, such as fpsetprec ( ) on FreeBSD, and 
f env on Linux. 
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Results of floating point arithmetics are stored in variables that may not offer the 
same level of precision the FPU is running at and in fact it is possible that none 
of the data types available matches a particular level of precision set on the FPU. 
Crucially, the precision setting of the FPU normally affects all floating point opera- 
tions on all floating point variables, regardless of type, e.g. information is lost when 
results are calculated with extended precision and stored in variables offering only 
single precision. A notorious error observed on systems which default to extended 
precision, in particular Linux on x86, occurs when comparisons between variables 
produce different outcomes depending on the position in the code — at one point 
the result calculated may still reside on the FPU and thus offer extended precision, 
whereas at a later point the result is truncated after being written to memory. This 
can lead to serious inconsistencies when data is held in an ordered tree. Compiler 
switches like -f float-store for gcc help in these cases. 

The commonly used gcc compiler offers three basic floating point types, float, 
double and long double, matching the three levels of precision mentioned above. 
The very nature of SOC means that observables span very many order of magni- 
tudes. If variables that accumulate results, such as moments, are too small (i.e. have 
a mantissa that is too small), smaller events may not accumulate at all any more 
once the variable has reached a sufficiently large value. This can skew estimates 
considerably where very large events occur very rarely. The macros flt_epsilon, 
dbl_epsilon and ldbl_epsilon in f loat.h give a suggestion of the relative scale 
of the problem. It can be mitigated by frequently "flushing" accumulating variables 
(see Sec. 1.3.5). 



1.3.4 Random Number Generators 

Random Number Generators (RNGs) are a key ingredient in many areas of compu- 
tational physics, in particular in Monte-Carlo and Molecular Dynamics simulations. 
The vast majority of them, strictly, are not random, but follow instead a determin- 
istic but convoluted computational path. RNGs are constantly being improved and 
evaluated, not least because of their use in cryptography. An introduction into the 
features of a good RNG can be found in the well-known Numerical Recipes (Press 
et al., 2007), with further details to be found in the review by Gentle (1998). 

A "good" random number generator is one that offers a reasonable compromise 
between two opposing demands, namely that of speed and that of quality. In most 
stochastic SOC models, the RNG is used very often and thus typically consumes 
about half of the overall CPU time. Improving the RNG is thus a particularly sim- 
ple way of improving the performance of an implementation. Because the variance 
(square of the standard deviation) of an estimate vanishes inversely proportional 
with the sample size it is based on, the performance of an implementation is best 
measured as the product of variance and CPU time spent "for it". However, one is 
ill-advised to cut corners by using a very fast RNG which has statistical flaws. The 



1 SOC computer simulations 



43 



resulting problem may be very subtle and might not show until after a very detailed 
analysis. 

One of the problems is the period of an RNG. Because RNGs generally have a 
finite state, they are bound to repeat a sequence of random numbers after a sufficient 
number of calls, at which point the simulation using the random numbers produces 
only copies of previous results. With improving hardware the RNG must therefore 
be re-assessed. A "good RNG" is a function of time, and very much a function of 
perception, as a mediocre RNG might appear to be a fantastic improvement over a 
poor RNG. It is good practise to use more than one random number generator to 
derive the same estimates and compare the results. 

The C library's implementation of rand ( ) is legendary for being unreliable and 
can be very poor. At the very least, it is essentially uncontrolled, although, of course, 
standards exist, which are, however, not always adhered to. It is fair to say that 
pure linear congruential RNGs are somewhat (out-)dated and indeed rarely used. 
They are, however, sometimes combined or enhanced with more sophisticated tech- 
niques. In recent years, the Mersenne Twister (Matsumoto and Nishimura, 1998; 
Matsumoto, 2008) has become very widely used, yet, criticised by Marsaglia (2005) 
who proposed in turn KISS (Marsaglia, 1999, but see Rose, 2011), which is a re- 
markably simple RNG. The GNU Scientific Library (Galassi et al, 2009) contains 
an excellent collection of random number generators. 

Somewhat more specific to the use of RNGs in SOC models is the frequent de- 
mand for random bits, for example in order to decide about the direction a particle is 
taking. Because every acceptable RNG is made up of equally random bits, each and 
everyone of them should be used for random booleans. These bits can be extracted 
one-by-one, by bit-shifting the random integer or by shifting a mask across, as in 

#define RNGJVIT JITS (32) 
#define RNG.TYPE unsigned long 

RNG.TYPE mt_bool_rand=0UL; 

RNG.TYPE mt_bool_mask=lUL< <(RNG_MT_BITS— 1); 

#define RNGJVIT .BOOLEAN ( ( mt.booljnask==(lUL<<(RNGJVIT.BITS- 1)) ) ? (( 
U mt_bool_mask=lUL, mt_bool_rand=genrancLint32()) & mt_bool_mask) : ( 
U mt_bool_rand & (mt_bool_mask+=mt_booLmask)) ) 

based on the Mersenne Twister. In general, bit shifts to the left using a+=a instead 
of a<<=i are faster, because the latter requires one more CPU cycle to write the 
constant l into the CPU's register. 

More generally, integer random numbers have to be chosen uniformly from the 
range {0, l,...,n— 1} suggesting the use of the modulo operation, r=rand()%n. 
However, if rand ( ) produces random integers uniformly from up to and including 
rand_max, then the modulo operation skews the frequencies with which random 
number occurs towards smaller values if rand_max+i is not an integer multiple of 
n. The effect is of order n/ (rand_max+i ) and thus is negligible if n is significantly 
smaller than rand_max. However, picking a site at random on a very large lattice 
or an element from a very large stack, this effects becomes a realistic concern. In 
that case, the modulo operation can be used on a random number drawn uniformly 



44 



Gunnar Pmessner 



among integers from up to and including R—l, where R is a multiple of n and 
ideally the largest multiple of n less or equal to rand_max+i: 

const long long int n=...; 

/* The constant multiple jninus _1 is made to have type as the return 
* value ofrandf). */ 

const int multiple_rainus_l=(n*((((long long int)RAND_MAX) + lLL)/n))- ILL; 
int r; 

#define RANDOM(a) while ((r=rand())>multiple_minus_l); (a)=r%n 

where muitipie_minus_i plays the role of R— 1. When determining the maxi- 
mum multiple, it is crucial that the operation rand_max+i is performed using a type 
where the addition does not lead to rounding or integer overflow. The latter is also 
the reason why one is subtracted in the expression for multiple_minus_l, which 
otherwise might not be representable in the same type as the return value of rand ( ) , 
which is necessary to avoid any unwanted type casting at run time. 16 

The initial seed of the RNG needs to be part of the output of the programme it 
is used in, so that the precise sequence of events can be reproduced in case an error 
occurs. Some authors suggest that the initial seed itself should be random, based, 
for example, on /dev/random, or the library functions time ( ) or clock ( ) , 17 and 
that the RNG carries out a "warm-up-cycle" of a few million calls (Jones, 2012). 
After that, it is sometimes argued, chances are that one sequence of (pseudo) ran- 
dom numbers is independent from another sequence of random numbers generated 
by the same RNG based on a different seed. Fortunately, some RNGs, in particular 
those designed for use on parallel machines, offer a facility to generate sequences 
that are guaranteed to be independent. Where poor-man's parallel computing (many 
instances of the same simulation running with different seeds) takes place, inde- 
pendent sequences are of much greater concern than in situations where different 
parameter settings are used in different instances. In the former case the data of all 
instances will be processed as a whole, probably under the assumption that it is ac- 
tually independent. In the latter case, the results will enter differently and using even 
an identical sequence of random numbers will probably not have a noticeable effect. 
All these caveats are put in perspective by the fact that most SOC models fed by a 
slightly differing sequences of pseudo random numbers take "very different turns in 
phase space" and thus will display very little correlations. 



16 This is one of the many good reasons to use constants rather than macros (van der Linden, 1994; 
Kernighan and Pike, 2002). 

17 Both functions are bad choices on clusters where several instances of the same programme are 
intended to run in parallel. The function time ( ) changes too slowly (returning the UN*X epoch 
time in seconds) and the function clock ( ) wraps after about 36 minutes, so that neither function 
guarantees unique seeds. In general, seeding is best done explicitly. 
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1.3.5 Output 

As mentioned above, it is generally advisable to output and flush data frequently 
in chunks, resetting accumulating variables afterwards. Even when output occurs 
every second, the overhead in terms of the CPU and real time spent by the system is 
likely to be negligibly small. 

Where data is written to a file in large quantities or frequently, buffered I/O as 
provided by stdio through the print f -family of library calls is usually much faster 
than writing immediately to the file using unistd's write. There are two caveats 
to this approach: Firstly, depending on the size of the buffer and thus the frequency 
of writing, a significant amount of CPU time may be lost if the program terminates 
unexpectedly. To avoid corrupt data, f flush ( ) should be used rather than allowing 
the buffer to empty whenever it reaches its high-water mark. Secondly, if buffering 
I/O has a significant impact on the computational performance, the data may better 
be processed on-the-fly rather than storing it in a file. 

In the following, stdio is used for its convenient formatting capabilities, pro- 
vided by the plethora of flags in the formatting string of a printf call. To avoid the 
problems mentioned above, buffers are either flushed after each chunk by means of 
f flush, or buffering is switched to buffering line by line, using setiinebuf . 

To avoid unexpected interference of the operating system with the simulation, 
operations should be avoided that can potentially fail because the environment 
changes. This applies, in particular, to read and write access to files. In any case, 
such operations need to be encapsulated in an if condition that catches failing sys- 
tem calls and triggers a suitable remedy. 

Output of chunks should therefore happen through the stdout stream which is by 
default open at the time of the program start. As the output is usually used in post- 
processing it needs to be retained, which can be achieved by re-directing stdout 
into a file. In the typical shell syntax this can be done in the command line by, say, 
./simulation > output.txt. To allow easy post-processing, every line should 
contain all relevant simulation parameters, such as the system size, the number of 
the chunk (a counter), the number of events per chunk, the initial seed of the random 
number generator (RNG), in fact, everything that is needed to reproduce that line 
from scratch or to plot the relevant (derived) data. Typical examples are moments to 
be plotted against the system size and moment ratios, involving different moments 
of the same observable. Using post-processing tools to wade through vast amounts 
of data to find the missing piece of information to amend a line of data can require 
significant effort and is highly error-prone. 

Repeating the same output (system size, RNG seed etc) over and over seemingly 
goes against the ethos of avoiding redundant information, which should be applied 
when setting up a computer simulation (to avoid clashes), but is wholly misplaced 
when it comes to data output. In fact, redundancy in output is a means to measure 
consistency and a matter of practicality as almost all basic post-processing tools are 
line-oriented. 

In some rare cases, an action by the simulation or an event on the system can 
result in a signal being sent to the running instance of the program. In response 
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the program suspends the current operation, executes a signal handler and contin- 
ues where it left off. In principle, the signal should not lead to inconsistent data or 
behaviour; in fact, it is probably the most basic but also a very convenient way to 
communicate with a running program. For example 

#include <signal.h> 

void sighup_handler(int signo); 

signal(SIGHUP, sighupjiandler); 

void sighup_handler(int signo) 
{ 

finish_asap=l; 
} 

assigns the signal handler sighup_handier to deal with the signal sighup, which 
can be sent to the program using kill -hup. 

There is a rare situation when the signal interrupts in a way that it leads to un- 
expected behaviour, namely when it arrives while a "slow system call" is executed, 
i.e. an operation that is performed by the kernel on behalf of the programme, but 
which can take a long time to complete, such as pause, sleep, but also write to 
so-called pipes. Without discussing the technical details of the latter, it can lead to 
inconsistencies in the output which might not be detected in the post-processing. For 
example, a chunk may contain truncated lines and thus may lack certain information 
or data, which the post-processing tools might treat as zeroes. Apart from a graphi- 
cal inspection of the data, two measures may therefore be advisable: Firstly, output 
can be encapsulated in calls of sigprocmask which allows temporary suspension 
of the delivery of signals. Secondly, a chunk can be terminated by a single line con- 
taining a keyword to indicate the successful completion of the output (i.e. without 
catching an error, in particular not an "interrupted system call", eintr), such as the 
tag (see below) #compieted. Simply counting the number of occurrences of that 
tag and comparing to (supposed) the number of valid chunks can pick up inconsis- 
tencies. In large scale simulations, where disk space can be a problem leading to 
truncated files as the system runs out of file space, this is particularly advisable. 

After a chunk has been written out, variables collecting data have to be reset. 
Where PDFs are estimated, sweeping across the entire histogram can become ex- 
pensive and therefore performing all relevant steps simultaneously is advantageous 
for the overall performance. Using one of the examples above (Sec. 1.2.1.1): 

long long total=0; 

for (i=0; i < SMALL2MEDIUM_THRESHOLD ; i++) 
if (histo_small[i]) { 
printf(...); 

total+=histo_small [i] ; 
histo_small[i]=0; 

} 

printf(''out_of_range:„%i\n", histo_out_of_range); 
total+=histo_out_of_range ; 



1 SOC computer simulations 



47 



histo_out_of_range=0; 
printf("total:„%lli", total); 

The final line allows the user to compare the number of events collected in the 
histogram to the number of events expected. It is a computationally cheap additional 
check for data consistency. 

To distinguish different types of output, such as moments of different observ- 
ables, data should be tagged by short keys that are easily filtered out in post pro- 
cessing. For example, if every line containing moments of avalanche sizes is tagged 
by #m_s i ze at the beginning, all relevant lines can be extracted very easily for exam- 
ple using grep ' ~#m_size' output . txt. To strip off the tags, one either appends 
I sed ' s/#m_size//' or includes the functionality of grep in the sed command, 

sed -n 's/"#M_SIZE//p' output.txt > output.txtJVLSIZE 

storing all relevant lines in output .txt_M_siZE for further processing by other 
tools. One very simple, but particularly powerful one is awk. For example, the av- 
erage across the seventh column starting with the 101st chunk (stored in the first 
column) can be calculated using 

awk ' { if ($1>100) {m0++; ml+=$7;} } END { printf ("%i %10.20g\n", mO, ml/mO); 
U } ' output.txtJVLSIZE 

All of this is very easily automated using powerful scripting languages (in particular 
shell scripts, awk, sed and grep), and more powerful (interpreted) programming 
languages, such as perl or python, which provide easy access to line-oriented data. 
In recent years, XML has become more popular to store simulation parameters as 
well as simulation results. 



1.4 Summary and conclusion 

The early life of SOC was all about computer models that showed the desired fea- 
tures of SOC: Intermittent behaviour (slow drive, fast relaxation) displaying scale 
invariance as observed in traditional critical phenomena without the need to tune 
a control parameter to a critical value. After many authors had (mostly with lit- 
tle success) attempted to populate the universality class of the BTW Sandpile, a 
range of SOC models was proposed firstly as a paradigm of alternative universal- 
ity classes and later to highlight specific aspects of SOC, such as non-conservation 
(as for example in the Forest-Fire Model), non-Abelianness (as for example in the 
Olami-Feder-Christensen Model) and stochasticity (as for example in the Manna 
Model). 

Many of these models have been studied extensively, accumulating hundreds of 
thousands of hours of CPU time in large-scale Monte Carlo simulations. A finite size 
scaling analysis of the data generally produces a set of two to eight exponents, which 
are supposedly universal. It turns out, however, that very few models display clean, 
robust scaling behaviour in the event size distribution, although it is remarkably 
broad for many models. 
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Of the models discussed above, the Manna Model displays the clearest signs 
of scale invariance. There is wide consensus that it is the same universality class 
as the Oslo Model (Christensen et al, 1996; Nakanishi and Sneppen, 1997). In the 
conservative limit and in the near-conservative regime, the Olami-Feder-Christensen 
Model also displays convincing moment scaling, but less so for smaller values of the 
level of conservation. Numerical artefacts may play a significant role in its scaling 
(Drossel, 2002). 

The Forest Fire Models is widely acknowledged for failing to display finite size 
scaling in the event size distribution (Grassberger, 2002; Pmessner and Jensen, 
2002a), although its moments still display some scaling (Pmessner and Jensen, 
2004). The contrast is even sharper in the Bak-Tang-Wiesenfeld Model: Some scal- 
ing is known analytically (Majumdar and Dhar, 1992; Ivashkevich, 1994; Ivashke- 
vich et al, 1994; Dhar and Manna, 1994), yet the event size distribution seems at 
best be governed by multiscaling (Tebaldi et al, 1999; Drossel, 1999, 2000; Dorn 
efa/.,2001) 

While analytical approaches receive increasing attention, numerical techniques 
remain indispensable in the development and analysis of models which are tailor- 
made to display specific features or to mimic experimental systems. Models de- 
veloped more recently are usually implemented in C, producing numerical data in 
Monte-Carlo simulations. It is fair to say that the careful data analysis requires as 
much attention to detail as the implementation of the model in the first place. 

While the classic data-collapse and more immediate tests for scaling dominated 
the early literature of SOC, more recently the finite size scaling of moments (Tebaldi 
et al, 1999) has become the predominant technique for the extraction of scaling 
exponents. Apart from identifying the mechanism of SOC, the main purpose of the 
numerical work is to establish universality and universality classes among models, 
as well as their relation to natural phenomena. One may hope that these efforts will 
eventually help to uncover the necessary and sufficient conditions for SOC. 

Acknowledgements The author gratefully acknowledges the kind support by EPSRC Mathemat- 
ics Platform grant EP/I0191 1 1/1. 



Appendix: Implementation details for binning 

To implement binning in computer simulations of SOC models it is advisable to per- 
form simple bit manipulations on basic, integer-valued observables. It often suffices 
to implemented three levels of coarse graining or less, for example 

#define SMALL2MEDIUM.THRESHOLD (1LL<<15) 
long long histo_small[SMALL2MEDIUM_THRESHOLD]={0LL}; 
#define MEDIUM2LARGE.THRESHOLD (1LL<<30) 
#define MEDIUM.SHIFT (12) 

long histo_medium[(MEDIUM2LARGE_THRESHOLD— 

U SMALL2MEDIUM_THRESHOLD)>>MEDIUM_SHIFT]={0L}; 

#define LARGE.THRESHOLD (1LL<<45) 
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#define LARGE.SHIFT (27) 

int histo Jarge [(LARGE.THRESHOLD— MEDIUM2LARGE_THRESHOLD) > > 
U LARGE_SHIFT]={0}; 
int histo_out_of_range=0; 
long long int s; /* event size */ 



if (s<SMALL2MEDIUM_THRESHOLD) histo.small[s]++; 

else if (s<MEDIUM2LARGE_THRESHOLD) histo_medium[(s- 

U SMALL2MEDIUM.THRESHOLD)> >MEDIUM.SHIFT]++; 

else if (s <LARGE_THRESHOLD) histo Jarge[(s—MEDIUM2LARGE_THRESHOLD) 

U >>LARGE.SHIFT]++; 

else histo_out_of_range++; 

Here the event size to be tallied is s. In the block of if statements, it is compared 
to various thresholds before it is rescaled and counted into a histogram. Because 
vectors in many programming languages start with index 0, a shift an offset is sub- 
tracted as well. It can pay of to re-arrange the if statements as to test against the 
most frequent case as early as possible. One case, in the present example the last 
one, counts the number of times the counter overspills, here histo_out_of_range. 

Some subtleties of the above implementation are worth discussing. Firstly, the 
types used for the histogram typically decrease in size with increasing event size 
while the size of the type needed to represent the event size at the respective thresh- 
olds increases. This is because normally the frequency is an inverse power law of 
the event size. Great care must be taken to avoid unnecessary typecasts and unde- 
sired outcomes, as some languages, in particular C, are rather idiosyncratic when it 
comes to (integer) type-promotion in comparisons, in particular when they involve 
signs. 

In the above examples, automatic vector variables are used and initialised by as- 
signing { } , which is expanded by the compiler to a suitable size by adding zeroes. 
Initialisation of vectors in C has been further simplified in the C99 standard. 

Secondly, it is important to choose the thresholds together with the planned bit- 
shifts, in order to avoid an off-by-one error. The problem is that, say, 

s <MED IUM2 LARGE_THRE SHOLD , does not imply 

(s— SMALL2MEDIUM_THRESHOLD)/((l<<MEDIUM_SHIFT) < ( 

U MEDIUM2LARGE_THRESHOLD— SMALL2MEDIUM_THRESHOLD)/( 1 < < 

U MEDIUM.SHIFT) 

because for some s<medium2large_threshold their bitshifted value 
s»medium_shift in fact equals medium2large_threshold»medium_shift, 
namely precisely when medium2large_threshold is not an integer multiple of 
k<medium_shift. It is therefore a matter of defencive programming to write the 
thresholds for the macros in this form: 

#defme MEDIUM2LARGE.THRESHOLD ((1LL«18) * (1LL<<MEDIUM_SHIFT)) 

As for a rudimentary output routine 

for (i=0; i < SM ALL2MEDIUM_THRESHOLD ; i++) 
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if (histo_small[i]) printf("%i„%i„%lli„%i\n", i, i, histo_small[i], 1); 
for (i=0; i<((MEDIUM2LARGE_THRESHOLD— SMALL2MEDIUM_THRESHOLD) 
U >>MEDIUM.SHIFT); i++) 

if (histo_medium[i]) printf("%li„%i„%li„%i\n", ((long) 
U SMALL2MEDIUM.THRESHOLD)+(((long)(i))<<MEDIUM.SHIFT), i, 
U histo_medium[i], 1 < <MEDIUM_SHIFT); 

for (i=0; i<((LARGE_THRESHOLD— MEDIUM2LARGE_THRESHOLD)>> 
U LARGE.SHIFT); i++) 

if (histo_large[i]) printf("%lli J&i„%i„%i\n", ((long long) 
U MEDIUM2LARGE_THRESHOLD)+(((long long)(i))< <LARGE.SHIFT), i, 
L, histojargep], 1 < <LARGE_SHIFT); 
printf("out_of_range:„%i\n", histo_out_of_range); 

care must again be taken that the formatting of the output is in line with the type 
of the data and does not spoil it. Fortunately, most modern compilers spot clashes 
between the formatting string used in print f and the actual argument. As discussed 
below, it is generally advisable to have only one output stream, namely stdout, 
and to use tags to mark up data for easy fetching by post-processing tools. In the 
example above, the bins have not been rescaled by their size which instead has been 
included explicitly in the output. A sample of the PDF can be inspected by plotting 
the third column divided by the fourth against the first. 
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